An efficient numerical model is developed in this work for the bifurcation analysis within Föppl-von Karman plate theory. This procedure is based on the use of the method of fundamental solutions in a high order continuation method that relies on Taylor series expansion with respect to a path parameter. This numerical model has an adaptive step length, which is effective especially for solving nonlinear problems and detecting bifurcation points. Despite of the huge application field of nonlinear elasticity, the method of fundamental solution was very rarely applied in this field and never to a nonlinear elastic plate model. The governing equations are strongly formulated in terms of the two unknowns: the deflection and the stress function. The singular value decomposition regularization is used to overcome the difficulty of the resulting system ill-conditioning. The accuracy and efficiency of the numerical model are illustrated on buckling numerical examples.