<p>This paper presents a hybrid technique for solving the partial integro-differential equation (PIDE) governing option pricing, a model rigorously derived from financial principles such as replication portfolios and risk-neutral valuation. The proposed numerical method combines the Crank–Nicolson finite difference scheme with an operational matrices approach based on Morgan–Voyce polynomials. First, the model is discretized in time using the trapezoidal rule. At each time step, the solution is approximated as a truncated series of Morgan–Voyce polynomials. The resulting approximated model and boundary conditions are then collocated, transforming the complex PIDE into an iterative linear system. The required initial condition for this iterative system is provided by projecting the terminal payoff onto the polynomial basis. Furthermore, the stability and convergence of the method are rigorously investigated. The effectiveness and accuracy of the hybrid scheme are demonstrated through five numerical test problems. The first three are manufactured solutions used to verify the method’s correctness, while the final two are real-world cases with no known analytical solutions, highlighting the scheme’s practical utility.</p>