We develop a general numerical method for the study of bifurcation problems (turning points, simple bifurcations, pitchfork bifurcations and a host of other singularities) that occur in the solutions to boundary value problems that depend on parameters. Consider time-evolution equations of the form [Formula: see text] where u∈Rn is the state vector, λ∈Rp is a vector of parameters and F is a sufficiently smooth function with values in Rn. If this system describes a time-dependent boundary value problem then n≫p and the Jacobian Fu typically has a special sparsity pattern (banded, block banded,…). A bordered matrix is an (n+m, n+m) matrix M (n≫m) that consists of a main (n, n) block A in the upper left corner (typically A is Fu or closely related to it) and bordering rows and columns. We consider methods for computing singularities of the steady state solutions F(u, λ)=0 in the general sense of Golubitsky & Schaeffer [1985]. In the case of the cusp catastrophe, Griewank & Reddien [1989] showed that Newton schemes can be set up to compute such singularities by solving only a small number of bordered linear systems in each Newton step. Their method extends similar, older methods for the simplest catastrophe, the turning point. We show that it is possible to set up a general scheme to compute any singularity with a small number of solutions of bordered systems in each Newton step; this number depends on the type of singularity but not on n.