The paper uses a multigrid technique for the numerical computation of bifurcation points of large sparse systems of nonlinear, parameterized equations arising in the discretization of certain elliptic boundary value problems. More specifically, problems of the form (1) \(Lu=f(\lambda,u)\) in \(\Omega\), \(Bu=0\) on \(\partial \Omega\) are considered where \(\Omega \subset R^ n\) is a bounded domain, L a linear elliptic differential operator, B the boundary condition with \(B0=0\), and \(f: R^ 1\times R^ 1\to R^ 1\) a smooth nonlinearity with \(f(\lambda,0)=0\). The approach is based on the classical Ljapunov-Schmidt procedure for simple bifurcation points on the trivial branch \(u=0\). The multigrid technique is applied to the iterative solution of the discretized form of the resulting linear equation \(L_ 0v=g\), \(v\in Z\), \(g\in R(L_ 0)\), where \(L_ 0=L-f_ u(\lambda,0)\), dim N(L\({}_ 0)=1\), codim R(L\({}_ 0)=1\), and Z is a complement of \(N(L_ 0)\). The method follows that proposed by \textit{W. Hackbusch} [ibid. 16, 201-215 (1979; Zbl 0403.65043)] for singular problems. It is combined with a nested iteration technique. The paper shows how this bifurcation iteration can be implemented in the case of the large sparse systems arising in the discretization of (1). Several numerical examples show that the method can be very efficient for various problems in \(R^ 1\) and \(R^ 2\).