I have transformed the BSM PDE $$\frac{\partial V}{\partial t} + \frac{\sigma^2}{2}S^2 \frac{\partial^2 V}{\partial S^2} + rS \frac{\partial V}{\partial S} - rV = 0 $$ to $u(\tau,x) = V(T-\tau,S_{0} e^{x})$ with a change of variables $x = \ln(S/S_{0})$ and $\tau = T -t$ to $$ \frac{\partial u}{\partial \tau} = \frac{\sigma^2}{2} \frac{\partial^2 u}{\partial x^2} + (r - \frac{\sigma^2}{2})\frac{\partial u}{\partial x} - ru $$ For the discretization I use backward difference for the time derivative and central differences for the spatial derivatives. So I have an implicit finite difference scheme of the form $$ A u_{i}^{n+1} = u_{i}^{n} $$ where $i$ is the index for the spatial grid (from 1 to $M = 2B/\Delta x$) and $n$ for the time grid (from 1 to $N = T/\Delta\tau$) and $$ A = \begin{pmatrix} \beta_{1} & \gamma_{1} & 0 & 0 & \cdots & 0\\ \alpha_{2} & \beta_{2} & \gamma_{2} & 0 & \cdots & 0 \\ 0 & \alpha_{3} & \beta_{3} & \gamma_{3} & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots &\ddots &\vdots \\ 0 &\cdots & & & \alpha_M & \beta_{M} \end{pmatrix} $$ and $\beta_i = 1 - k (\frac{\sigma^2}{\Delta x^2} - r)$, $\gamma_i = k(- \frac{\sigma^2}{2 \Delta x^2} - (r - \frac{\sigma^2}{2}) \frac{1}{2\Delta x})$ and $\alpha_i = k( -\frac{\sigma^2}{2 \Delta x^2} + (r - \frac{\sigma^2}{2})\frac{1}{2 \Delta x})$ $\forall 2 \leq i \leq N-1$ From the boundary conditions for a European call $$\lim_{S \rightarrow 0} V(S,t) = 0 \qquad \lim_{x \rightarrow -\infty} u(x,\tau) = 0$$ I set $\beta_1 = 1$, $\gamma_{1} = 0$ and $u_{i}^{n} = 0$ as the lower bound and from $$\lim_{S \rightarrow \infty} V(S,t) = S - K e^{-r(T-t)} \qquad \lim_{x \rightarrow \infty} u(x,\tau) = S_{0} e^{x} - K e^{-r(T-t)}$$ I set $\alpha_M = 0$, $\beta_{M} = 1$ and $u_{M}^{n} = S_{0} e^{x_{M}} - K e^{-r(T-t)}$ to solve this scheme.
Now I want to use von Neumann boundaries and I'm not sure how it's done in the transformed PDE. So far I have tried it by leaving the lower boundary as it is and for the upper boundary I used $$\lim_{S \rightarrow \infty} \frac{\partial V(S,t)}{\partial S} = 1 \qquad \lim_{x \rightarrow \infty} \frac{\partial u(x,\tau)}{\partial x} = S_{0} e^{x} $$ From the central difference of the first spatial derivative I get then $$ \frac{u_{M+1}^{n+1} - u_{M-1}^{n+1}}{2 \Delta x} = e^{x_{M}} S_{0} \\ u_{M+1}^{n+1} = 2 \Delta x e^{x_{M}} S_{0} - u_{M-1}^{n+1}$$ Now I insert this into the last row of the scheme $$\alpha_{M} u_{M-1}^{n+1} + \beta_{M} u_{M} + \gamma_{M} (2 \Delta x e^{x_{M}} S_{0} - u_{M-1}^{n+1}) = u_{M}^{n} \\ (\alpha_{M} - \gamma_{M}) u_{M-1}^{n+1} + \beta_{M} u_{M} = u_{M}^{n} - 2 \Delta x e^{x_{M}} S_{0} \gamma_{M} $$ So the boundary in the scheme becomes $\alpha_{M} = (\alpha_{M} - \gamma_{M})$, $\beta_M = \beta_M$ and $u_{M}^{n} = u_{M}^{n} - 2 \Delta x e^{x_{M}} S_{0} \gamma_{M}$
When I implement this in Matlab with the Dirichlet boundary I get good results, but something must be wrong with the von Neumann implementation because the results for large $S$ have very large errors. Can someone please help me?
Thanks a lot