We adapt numerical continuation methods to compute all solutions of finite difference discretizations of nonlinear boundary value problems involving the Laplacian in two dimensions. New solutions on finer meshes are obtained from solutions on coarser meshes using a complex homotopy deformation. Two difficulties arise. First, the number of solutions typically grows with the number of mesh points and some form of filtering becomes necessary. Secondly, bifurcations may occur along homotopy paths of solutions and efficient methods to swap branches are developed when the mappings are analytic. For polynomial nonlinearities we generalize an earlier strategy for finding all solutions of two-point boundary value problems in one dimension and then introduce exclusion algorithms to extend the method to general nonlinearities.