There is an example of using my non-smooth suggestion with NDSolve (including code) in my Vol II book, pages 462-464. It works quite well, as Fig 10.5 there shows.

This looks very good, Alan.

I reckom NDSolve (using Bulisch-Stoer) can handle many kinds of pesky deltas:

1. avoiding oscillations and negative values

2. what is accuracy order (high?), order 4 in space, order 5 in time?

3. sanity check: integral == 1. (cdf)

4. what about small and large T?

Now, I have modelled the PDE for option gamma (differentiate BSPDE twice wrt S) with a Gaussian as delta. I used C++, centred differencing in S and Bulisch-Stoer in time (Boost C++ odeint library). I get excellent results.

Will post details later today.

Meanwhile, back at the ranch, people are looking for fixes for traditional fdm. It is based on flawed assumptions and expectations. A kind of wasted effort.

They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.