In this paper we present a new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation. The scheme uses a symmetrical multi-point difference formula to represent the partial differentials of the two-dimensional variables, which can improve the accuracy of the numerical solutions to the order of Δ x 2 N q + 2 when a ( 2 N q + 1 ) -point formula is used for any positive integer N q with Δ x = Δ y , while N q = 1 equivalent to the traditional scheme. On the other hand, the new scheme keeps the same form of the traditional matrix equation so that the standard algebraic eigenvalue algorithm with a real, symmetric, large sparse matrix is still applicable. Therefore, for the same dimension, only a little more CPU time than the traditional one should be used for diagonalizing the matrix. The numerical examples of the two-dimensional harmonic oscillator and the two-dimensional Henon–Heiles potential demonstrate that by using the new method, the error in the numerical solutions can be reduced steadily and extensively through the increase of N q , which is more efficient than the traditional methods through the decrease of the step size.