<p>This study presents a novel and accurate numerical scheme for solving the two-step option pricing partial differential equation involving two spatial and one temporal independent variables. The complexity of the model, arising from the multidimensional nature of financial instruments, necessitates the development of efficient and precise computational techniques. To this end, we introduce an innovative approximation method based on a quasi-interpolating cubic spline for bivariate functions, which offers superior smoothness and flexibility in handling the spatial components of the problem domain. For the temporal approximation, we employ a parametric finite difference scheme designed to enhance both accuracy and stability in capturing the system’s dynamics. The proposed method is rigorously analyzed from a theoretical standpoint, with detailed proofs ensuring convergence, thereby establishing a solid mathematical foundation. Extensive numerical experiments are conducted on representative test cases to evaluate the performance of the proposed scheme. The results confirm that our method achieves high accuracy and demonstrates its robustness and effectiveness in solving multidimensional models in financial applications</p>