Many mathematical models in biology and physiology are represented by systems of nonlinear differential equations. In recent years these models have become increasingly large-scale and multiphysics, as increasing amounts of data are available on the properties and behaviour of biological systems. Often an observed behaviour of interest in a model may be written as a linear functional. A key question therefore is to determine which terms in the model have the greatest effect on functionals of interest. An approach for answering this question has recently been developed, called model reduction using a posteriori analysis. The method was initially developed for systems of nonlinear initial value ordinary differential equations, and automatically identifies regions of the computational domain and components of the model solution where an accurate mathematical representation of the model is required to accurately calculate a linear functional of interest. Initial-value ordinary differential equations can be written as a first-order derivative term plus an algebraic 'reaction' term. In previous work on model reduction using a posteriori analysis the algebraic 'reaction' term is removed from the equations in the reduced model. The first contribution of this thesis is to extend the method so that the first-order derivative term is removed from the differential equation instead of the algebraic 'reaction' term, resulting in a quasi-steady state approximation in automatically identified regions of the domain and components of the solution. The second contribution of this thesis is to extend the method to boundary value problems and partial differential equations. We consider differential equations with algebraic terms, first order terms and second order terms, any combination of which may be nonlinear. The method is used to automatically simplify several differential equation models including models of chemotaxis and tissue-level cardiac electrophysiology.