Abstract.Computation iteration schemes and memory allocation technique for finite difference method were presented in this paper. The transformed form of a groundwater flow problem in the generalized curvilinear coordinates was taken to be the illustrating example and a 3-dimensional second order accurate 19-point scheme was presented. Traditional element-byelement methods (e.g. SOR) are preferred since it is simple and memory efficient but time consuming in computation. For efficient memory allocation, an index method was presented to store the sparse non-symmetric matrix of the problem. For computations, conjugate-gradientlike methods were reported to be computationally efficient. Among them, using incomplete Choleski decomposition as preconditioner was reported to be good method for iteration convergence. In general, the developed index method in this paper has the following advantages: (1) adaptable to various governing and boundary conditions, (2) flexible for higher order approximation, (3) independence of problem dimension, (4) efficient for complex problems when global matrix is not symmetric, (5) convenience for general sparse matrices, (6) computationally efficient in the most time consuming procedure of matrix multiplication, and (7) applicable to any developed matrix solver.