I don't like these BCs and I suspect the near-field BC also not correct. Fichera at S= 0 gives du/dt + ru = 0.A trick that I was using for this kind of problem: try using the boundary conditions [$]u^{N} = 2 u^{N-1} - u^{N-2} [$], [$]u^{0} = 2 u^{1} - u^{2}[$] modelling [$]\partial_{xx}u=0[$] (linear behavior at infinity). AFAIR, this saved me to compute complex, payoff dependent, Dirichlet or Neuman boundary conditions. Obviously it only works for linear at infinity type payoff.This is the plot with spot prices at time T (vector ST in the code) on x-axis and option prices at time 0 ((F-K)*exp(-r*T) in the code) on y-axis. As you can see the there is a problem when the option price approaches the biggest value of the spot price, since the curve goes from linear to exponential, I guess this is due to the fact that the right end condition is missing
It's mostly old hat at this stage.