Numerical delta of Bond Options

I'm trying to calculate the delta for bond Call options. I'm using the vasicek model which gives the following solution for a Zero-coupon bond call option:

$Z = N P(t,S) \Phi(d_1) - K P(t,T) \Phi(d_2).$

Differentiating the above equation with respect to the underlying $P(t,S)$ gives the classical delta

$\Delta = N \Phi(d_1).$

where

$d_1 = \frac{ln \left(\frac{NP(t,S)}{KP(t,T)}\right) + \frac{\sigma_p^2}{2}}{\sigma_p}$ , $d_2 = d_1 - \sigma_p$.

Comparing the analytical delta with the numerical first derivative (with the call formula as input) using a central difference scheme

$Z' = \frac{Z(r_{i+1})-Z(r_{i-1})}{P(t,S)(r_{i+1})-P(t,S)(r_{i-1})}$

or any other higher-order stencil gives a very different result. For example, while the former approaches 1 as the price of P(t,S) increases the latter doesn't reach 0.7. See the figure below. The X-axis is the price of the bond $P(t,S)$. Here $t<T<S$.

These discrepancies increase whith the Option maturity.

Is there any trick to find numerical deltas for IR derivatives? I faced the same problem with Swaptions.

I think my confusion lies in the underlying.

Thank you. Allan

$\Delta = \frac{ \frac{\partial Z}{\partial r}-\frac{\partial P(t,T)}{\partial r}\frac{Z}{P(t,T)} } {\frac{\partial P(t,S)}{\partial r}}$