A rather general question regarding a specific problem I am facing with my Matlab implementation of the implicit FD method for this PDE:
\begin{equation} \frac{\sigma_s^2}{2}\frac{\partial^2 V}{\partial S^2} + \rho(t) \sigma_S \sigma_\alpha\frac{\partial^2 V}{\partial S \partial \alpha} + \frac{\sigma_\alpha^2}{2}\frac{\partial^2 V}{\partial \alpha^2} + \mu_s \frac{\partial V}{\partial S} + \mu_\alpha \frac{\partial V}{\partial \alpha} + \frac{\partial V}{\partial t} - r(t)V \end{equation}
When running my code, I obtain negative and non-monotonic-increasing values. However, they should be all positive and monotonic-increasing throughout the two-dimensions of the grid ($F$ and $\alpha$). I checked my boundary conditions again and again but they seem fine to me. Same holds for FD discretization of the PDE.
My question is, what could be the cause of my problem?
For example,
Is there a possibility that FD formula and boundary conditions are correct, but that a particular choice of step sizes is ruining the method?
How can I know if the problem is coding-related or that, for example, my FD matrix (A) is misspecified?
Am I correct that FD matrix A is not dependent on time, hence can be calculated outside the FD loop?
A more detailed version of my question is given here
Edit June 7: Extended the code to handle larger grids. Applied variable transformation to both F and vol (alpha) similar to here formula (2.16). (Obviously, the upper bound on F/x changed accordingly to account for transformation from F to x). I now obtain 'nicer' results, in the sense that values are positive and for large part of the grid monotonic-increasing. However, when $\beta=1$ and $\nu=0$, I do not obtain similar results as BS. Also, for the lower part of the $V$ matrix (upper part F/x domain), values are no longer monotonic increasing (along $\alpha$/column-dimension), see left down corner V matrix (near $F_{max}$ and $\alpha_{min}$) and right down corner V matrix (near $F_{max}$ and $\alpha_{max}$). $A$ does seem to satisfy the stability constraint: $||A^{−1}||_{\infty} = 1 \leq 1$. I have noted that the diagonals of $A$ tend to increase rapidly suddenly (to an order of magnitude +04) and then decrease, after it increases again, etc (example). Is this how $A$ would normally behave (after variable transformation)?