Abstract: To search direction of each solution branch at pitch-fork bifurcation point, a numerical method is proposed in this paper. Because this method does not require information of derivatives of tangential stiffness matrix, it is suitable for general nonlinear finite element program. Stability of equilibrium corresponding to each solution curve is judged by Lagrange-Dirichlet theorem. The solution path in post-buckling of a square plate is traced globally and equilibrium stability of the plate was judged. The stable equilibrium can simulate the whole process of snap-through of square plate during loading and unloading. The result in this paper agrees with available experimental data. In part I, numerical methods and algorithm are proposed. In part II, numerical results of snap-through of square plate are compared with experimental data.