I wrote this up, triggered by Billy7's nice remark.some test output
/* S exact delta CN delta Spline delta exact gamma CN gamma Spline gamma
59 0.330935233 0.330702564 0.330478248 0.040967715 0.04098453 0.041068871
60 0.372482798 0.372232265 0.372093554 0.042043376 0.042074873 0.042161739
61 0.41484882 0.414593361 0.414541127 0.042603904 0.042647319 0.042733408
62 0.457519085 0.457270132 0.457302102 0.042654217 0.042706223 0.042788541
63 0.499993235 0.499760129 0.499871257 0.042216737 0.04227377 0.042349769
64 0.541801022 0.541590678 0.541773643 0.041328768 0.041387329 0.041455004
65 0.582515647 0.582332471 0.582578239 0.040039354 0.040096257 0.040154187
66 0.62176381 0.621609821 0.621908227 0.038405882 0.038458444 0.038505789
67 0.659232357 0.659107471 0.659447781 0.036490702 0.036536856 0.03657332
68 0.694671633 0.694574052 0.694945475 0.034357967 0.034396306 0.034422069
69 0.727895851 0.727822509 0.728214631 0.032070838 0.032100608 0.032116243
70 0.758780913 0.758727922 0.759131055 0.029689185 0.029710218 0.029716604
71 0.787260172 0.787223244 0.787628682 0.027267796 0.027280425 0.027278649*/
Table.1 Output from Steps 1-4; data is K = 65, T = 0.25, r = b = 0.8, v = 0.3,NS = 325, NT = 1000
//
Examining the accuracy of the results in Table 1 we see that both the Crank Nicolson and cubic spline interpolator produce similar results for delta while the interpolator’s accuracy for gamma is less accurate than the Crank Nicolson method. What is the reason? For smooth functions centred divided difference approximation to the first and second derivatives result in second order accuracy. This approach is based on Taylor expansions. Cubic splines are polynomials and they are a third-order accurate approximation to a smooth function. The first derivative of the spline is a second-order approximation to the exact delta while the second derivative of the spline is a first-order approximation to the exact gamma. This theoretical result is proven mathematically in Alberg, Nilson and Walsh 1967, Theorem 3.11.1, page 94. But cubic spline interpolators can overshoot in order to always create a visually pleasing curves and this might account for the fact that the values for gamma are larger than the exact gamma or even the Crank Nicolson gamma.
When approximations to the second derivatives of a smooth function are needed, we can first approximate it by using the values of the first derivative as input in order to compute an approximation to the second derivative. This is called spline-on-spline computation (Hildebrand 1974, page 492 and Alberg, Nilson and Walsh 1967, pages 48-50) and the improvements in accuracy are visible. In the current case a useful exercise would be to compare the values for gamma based on this variant and the original approach taken to produce the results in Table 1.