The order of accuracy for my Crank-Nicolson scheme is wrong.

I've been trying to figure out what the issue is in my numerical scheme I've written to solve the particle in the box problem for the Schrodinger equation. I've used a Crank-Nicolson scheme to implicitly solve for the wave function numerically, the issue I'm running into is that the Crank-Nicolson scheme is second order accurate in space, but I'm only finding that my numerical solution is on first order accurate. I am near 100% positive that the Crank-Nicolson scheme is order (2,2) for the Schrodinger equation, just like the heat equation, and that there should be no special exceptions here. Because of this, I'm assuming there's an error somewhere in my code for solving of the numerical solution.

My function for solving the particle in a box is here:

function [wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N)
    m = ceil((xR - xL)/h); % Number of x-axis intervals
    h = (xR - xL)/m;
    n = ceil(T/k); % Number of t-axis intervals
    k = T/n;
    mu = k/h^2;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % natural units ---> MASS = 1, hbar = 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    DiagTerm = -2*(1 - 1i*2/mu);
    
    A = diag(DiagTerm*ones(1,m-1)) + diag(ones(1,m-2),1) + diag(ones(1,m-2),-1);
    
    B = -conj(A);
    
    wave_function = zeros(n+1,m-1);
    
    x_val = xL + h*(1:m-1);
    wave_function(1,:) = sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (x_val - xL));
    
    for tstep = 2:(n+1)
        wave_function(tstep,:) = ((A'*A)\A'*B * wave_function(tstep-1,:)')';
    end
    
    wave_function = [zeros(n+1,1), wave_function, zeros(n+1,1)]; % <-- Concatenates Columns to the left and right of wave_function for boundary conditions
    
    wave_function_final = wave_function(end,:);
    
    En = (N*pi)^2/8;
    
    error = sqrt(sum( abs( exp(-1i*En*T)*wave_function(1,2:end-1) - wave_function_final(2:end-1) ).^2));
end

And my live script I use for running is here:

Ms = [5 10 20 40 80];

T = 16/pi;
xL = -1;
xR = 1;
N = 1;

orders = zeros(1,length(Ms)-3);
Errors = zeros(1,length(Ms));

xx = -1:0.05:1;
En = (N*pi)^2/8;
p = plot(xx,real(exp(-1i*En*T)*sqrt(2/(xR - xL))*sin(N*pi/(xR-xL) * (xx - xL)) ), 'Color', [0 0 0], 'LineStyle',"--", 'LineWidth',2);
p.Color(4) = 0.2;

for j = length(Ms):-1:1
    disp(j)
    h = 2/Ms(j);
    k = 1/Ms(j)^2;
    
    [wave_function_final, wave_function, error] = Schrodinger_CrankNicolson_PIAB(h,k,T,xL,xR,N);
    hold on
    xx = -1:h:1;
    plot(xx,real(wave_function_final),'LineWidth',1.5)
    Errors(j) = error;
end

for j = 4:length(Ms)
    orders(j-3) = log(abs( (Errors(j) - Errors(j-1))/(Errors(j-1) - Errors(j-2)) )) / log( abs( (Errors(j-1) - Errors(j-2))/(Errors(j-2) - Errors(j-3)) ) );
end

grid on
hold off

disp(Errors)
disp(orders)
figure 
plot(Errors)