Abstract:The Crank-Nicolson method can be used to solve the Black-Scholes partial differential equation in one-dimension when both accuracy and stability is of concern. In multi-dimensions, however, discretizing the computational grid with a Crank-Nicolson scheme requires significantly large storage compared to the widely adopted Operator Splitting Method (OSM). We found that symmetrizing the system of equations resulting from the Crank-Nicolson discretization help us to use the standard pre-conditioner for the iterative matrix solver and reduces the number of iterations to get an accurate option values. In addition, the number of iterations that is required to solve the preconditioned system, resulting from the proposed iterative Crank-Nicolson scheme, does not grow with the size of the system. Thus, we can effectively reduce the order of complexity in multidimensional option pricing. The numerical results are compared to the one with implicit Operator Splitting Method (OSM) to show the effectiveness.