The engine makes use of a finite difference approach to solve the generalized Black Scholes model. The first step is to reformulate the equation to evolve backwards in time such that the known final state (the payoff) can be used as the initial condition. To do this set
\[\tau = T-t\]
such that
\[\frac{\partial U}{\partial \tau} = - \frac{\partial U}{\partial t}\]
Which leads to
\[\frac{\partial U}{\partial \tau}=\frac{1}{2} \sigma^2(S,\tau) S^2 \frac{\partial^2 U}{\partial S^2} + r(\tau)S \frac{\partial U}{\partial S} - r(\tau)U\]
In order to facilitate the use of a standardize finite difference solver, the Black Scholes equation is reformulated into a standard form (similar to the Heat Equation) by making the substitution
\[x = log(S)\]
Following the standard variable substitution technique
\[\frac{\partial U}{\partial S} = \frac{\partial U}{\partial x} \frac{\partial x}{\partial S} = \frac{1}{S}\frac{\partial U}{\partial x}\]
and
\[\frac{\partial^2 U}{\partial S^2} = \frac{1}{S}\frac{\partial}{\partial x}\Bigg(\frac{1}{S}\frac{\partial U}{\partial x}\Bigg)\]
which can be evaulated using the product rule and knowing that
\[\frac{1}{S} = e^{-x}\]
To give
\[\frac{\partial^2 U}{\partial S^2} = \frac{1}{S^2}\frac{\partial^2 U}{\partial x^2} - \frac{1}{S^2}\frac{\partial U}{\partial x}\]
Plugging this into the backwards evolving equation (the S, t dependence is herein implied)
\[\frac{\partial U}{\partial \tau}=\frac{1}{2} \sigma^2 S^2 \Bigg (\frac{1}{S^2}\frac{\partial^2 U}{\partial x^2} - \frac{1}{S^2}\frac{\partial U}{\partial x} \Bigg) \frac{\partial^2 U}{\partial S^2} + rS \frac{\partial U}{\partial S} - rU\]
which simplifies to
\[\frac{\partial U}{\partial \tau}=\frac{1}{2} \sigma^2 \frac{\partial^2 U}{\partial x^2} + (r-\frac{1}{2} \sigma^2)\frac{\partial U}{\partial x} - rU\]
This is clearly of the form
\[\frac{\partial U}{\partial \tau}=A \frac{\partial^2 U}{\partial x^2} + B\frac{\partial U}{\partial x} +CU\]
with
\[A = \frac{1}{2} \sigma^2 , B= (r-\frac{1}{2} \sigma^2) , C=-r\]
This transformed equation is then approximated by standard central difference equations but allowing for non-uniform spatial discretization.
\[\frac{\partial U}{\partial x}=\frac{U_{i+1} - U_{i-1}}{h_i + h_{i+1}}\]
\[\frac{\partial^2 U}{\partial x^2}=\frac{1}{h_i + h_{i+1}}\Bigg(\frac{U_{i+1} - U_{i}}{h_{i+1}} - \frac{U_{i} - U_{i-1}}{h_{i}}\Bigg)\]
where
\[h_i = x_i - x_{i-1}\]
The tempororal discretization is simply
\[\frac{\partial U}{\partial \tau}=\frac{U^{j+1} - U^{j}}{\tau}\]
where the time-step is a linear grid from 0 to T (maturity date) in M steps.