Verification of order of convergence of Implicit Euler Method to numerically solve Black-Scholes PDE

I’m trying to verify the order of convergence for implicit Euler method to numerically solve Black-Scholes PDE. Theory says that it should be $O(Delta t + Delta S^2).$ My code is working absolutely fine. I’m getting almost zero error when I compare my solution to the one obtained by Black-Scholes formula. However, I’m getting incorrect order of convergence. Does it have something to do with how I choose parameters? Could someone have a look at my code and tell me where I’m going wrong if there is a mistake somewhere? Or perhaps convey to me what option parameters need to be chosen so that I get the correct order of convergence?

Here is the MATLAB code:

function [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS)
% [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS)
% Prices European put or call option using implicit finite difference method.
% -- Input arguments --
%  pc = 'put' or 'call'
%  K  = strike price
%  T  = expiry time
%  r  = risk-free interest rate
%  sigma = volatility
%  Smax = max value of S used in finite difference grid
%  Nt = number of time steps
%  NS = number of grid points along S axis
% -- Output arguments --
%  V  = values of option at asset values in S
%  t  = vector of time points
%  S  = vector of asset prices

% Time discretization
Dt = T / Nt;
t = linspace(0, T, Nt+1);

% Asset price discretization
DS = Smax / (NS);
S = linspace(0, Smax, NS+1)';

% Storage for option values
V = zeros(length(S), length(t));

% Set initial condition at expiry and boundary values
switch lower(pc(1))
case 'p' % put option
  V(1,:) = K * exp(-r*(T-t));
  V(end,:) = 0;
  V(1:NS-1,end) = max( K - S(2:end-1), 0);
case 'c' % call option
  V(1,:) = 0;
  V(end,:) = Smax - K * exp(-r*(T-t));
  V(2:end-1,end) = max( S(2:end-1) - K, 0);
otherwise
  error(['unknown option type ' pc])
end

% Define index vector of interior asset prices
J = 2:NS;

c1 = sigma^2 * S(J).^2 * Dt / ( 2 * DS^2 );
c2 = r * S(J) * Dt / ( 2 * DS );

alpha = -c1 + c2;
beta  = 1 + r*Dt + 2*c1;
gamma = -c1 - c2;

% Create and factor the sparse matrix for the linear system.
A = spdiags([[alpha(2:NS-1); 0], beta, [0;gamma(1:NS-2)]], [-1, 0, 1], NS-1, NS-1);
[L,U,p] = lu(A,'vector');

% Time step backwards from expiry
for n = Nt:-1:1
  % Set up RHS vector from values at time step n+1
  b = V(J,n+1);
  % Adjust for boundary values at S = 0 and S = Smax
  b(1) = b(1) - alpha(1)*V(1,n);
  b(NS-1) = b(NS-1) - gamma(NS-1)*V(NS+1,n);

  % Solve linear system A * V(J,n) = b using existing LU factorization
  % V(J,n) = A  b; would recalculate factorization at each time step
  V(J,n) = U  ( L  b(p) );
end

function [c, dcds] = blackscholes(S, K, r, sigma, Tmt)
% [c, dcds] = blackscholes(S, K, r, sigma, Tmt)
% Black and Scholes formula for the value of a call option
% and its derivative with respect to volatility sigma
% S   = underlying asset price
% K   = strike price
% r   = risk-free interest rate
% Tmt = time to maturity = T - t where T = expiry
% If sigma is a vector of volatilities, then both the
% call value and its derivatives are vectors of the same size.
%
% Uses normpdf and normcdf from Statistics toolbox.

s = sigma * sqrt(Tmt);

d1 = ( log(S/K) + ( r + sigma.^2/2)*(Tmt) ) ./ s;
d2 = d1 -  s;

% Use normpdf and normcdf from Statistics toolbox
c = S .* normcdf(d1) - K * exp(-r*Tmt) * normcdf(d2);

% Derivative of call value w.r.t. volatility sigma
dcds = S .* normpdf(d1) * sqrt(Tmt);

% Test script for the functions to solve Black - Scholes PDE
% for the value of a call or put option.
% Define option parameters
r = 0.1;
sigma = 0.4;
T = 5/12;
K = 50;
pc = 'call';

% Define discretization parameters
% Asset price S goes from 0 to Smax
Smax = 100;

% Compare with Black and Scholes formula for call
Tmt = T;
[c, dcds] = blackscholes(S, K, r, sigma, Tmt);
Err = c - V(:,1);

fprintf("nConvergence with respect to Dtn")

fprintf("n%6s %10s %10s %8snn","NS", "Nt", "max error", "rate")

nrows = 8;
NS = 100;
Nt = 200;
err = zeros(nrows,1);
for row = 1:nrows
    Nt = 2 * Nt;
    [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS);
    for n = 1:Nt+1
        err_n = norm( V(:,n) - blackscholes(S, K, r, sigma, T - t(n)), Inf);
        if err_n > err(row)
            err(row) = err_n;
        end
    end
    if row == 1
        fprintf("%6d %10d %10.2en", NS, Nt, err(row))
    else
        ratio = err(row-1) / err(row);
        rate = log2(ratio);
        fprintf("%6d %10d %10.2e %8.3fn", NS, Nt, err(row), rate)
    end
end

fprintf("nConvergence with respect to Dsn")

fprintf("n%6s %10s %10s %8snn","Ns", "Nt", "max error", "rate")

nrows = 5;
NS = 50;
Nt = 2000;
err = zeros(nrows,1);
for row = 1:nrows
    NS = 2 * NS;
    %[U, x, t] = iEuler(a, L, T, f, gamma0, gammaL, u0, NS, Nt);
    [V, t, S] = bspdeimp(pc, K, T, r, sigma, Smax, Nt, NS);
    for n = 1:Nt+1
        err_n = norm( V(:,n) - blackscholes(S, K, r, sigma, T - t(n)), Inf);
        if err_n > err(row)
            err(row) = err_n;
        end
    end
    if row == 1
        fprintf("%6d %10d %10.2en", NS, Nt, err(row))
    else
        ratio = err(row-1) / err(row);
        rate = log2(ratio);
        fprintf("%6d %10d %10.2e %8.3fn", NS, Nt, err(row), rate)
    end
end
```

    Convergence with respect to Dt

    NS         Nt  max error     rate

   100        400   2.59e-04
   100        800   1.19e-04    1.120
   100       1600   5.79e-05    1.042
   100       3200   2.83e-05    1.032
   100       6400   1.40e-05    1.017
   100      12800   6.95e-06    1.009
   100      25600   3.46e-06    1.005
   100      51200   1.73e-06    1.002

Convergence with respect to Ds

    Ns         Nt  max error     rate

   100       2000   9.17e-04
   200       2000   2.53e-04    1.858
   400       2000   7.03e-05    1.848
   800       2000   2.20e-05    1.678
  1600       2000   8.92e-06    1.300