Can't get proper numerical convergence for complicated Advection-​Diffusion-​Reaction PDE

I have written an advection-diffusion-reaction PDE using a Crank-Nicolson finite difference scheme. The detail of which and the derivation can be found here: My Stack Exchange Question. I'm writing this to see if anyone see's an issue with the code I've written to solve this problem. The functions I have are

function Ctilde = myCDR(X,H,m,n,xs,zs,Cs,vd,alpha_p,Kr,C0,Q,zanchor,Stilde)
     K = @(zvar) K_func(zvar,alpha_p,Kr);
    dK = @(zvar) dK_func(zvar,alpha_p,Kr);
     Qtilde = @(xvar) zs*Q(xs*xvar)/Cs;
    dKtilde = @(zvar) dK(zs*zvar+zanchor);
    u0tilde = @(xvar) zs^2/xs*u0_func(xs*xvar);
     Ktilde = @(zvar) sign(zs*zvar+zanchor).*K(abs(zs*zvar+zanchor));
    sDtilde = @(zvar) sD(zs*zvar+zanchor);
    sAtilde = @(zvar) zs*sA(zs*zvar+zanchor);
    sRtilde = @(zvar) zs^2*sR(zs*zvar+zanchor);
    sDbelow = @(zvar) sD_below(zs*zvar+zanchor);
    sAbelow = @(zvar) zs*sA_below(zs*zvar+zanchor);
    Htilde = (H-zanchor)/zs;
    dztilde = Htilde/n;
    ztilde = 0:dztilde:Htilde;
    Xtilde = X/xs;
    dxtilde = Xtilde/m;
    xtilde = 0:dxtilde:Xtilde;
    gammaTerm = zeros(1,n+1);
    betaTerm = zeros(1,n+1);
    alphaTerm = 2*sDtilde(ztilde+0.5*dztilde) ...
        + dztilde*sAtilde(ztilde+0.5*dztilde) ...
        + dztilde^2*sRtilde(ztilde);
    betaTerm(2:end) = dztilde*(sAtilde(ztilde(2:end)+0.5*dztilde) - sAtilde(ztilde(2:end)-0.5*dztilde)) ...
        - 2*(sDtilde(ztilde(2:end)+0.5*dztilde) + sDtilde(ztilde(2:end)-0.5*dztilde));
    betaTerm(1) = dztilde*(sAtilde(0.5*dztilde) - sAbelow(-0.5*dztilde)) ...
        - 2*(sDtilde(0.5*dztilde) + sDbelow(-0.5*dztilde));
    gammaTerm(2:end) = 2*sDtilde(ztilde(2:end)-0.5*dztilde) ...
        - dztilde*sAtilde(ztilde(2:end)-0.5*dztilde) ...
        + dztilde^2*sRtilde(ztilde(2:end));
    gammaTerm(1) = 2*sDbelow(-0.5*dztilde) ...
        - dztilde*sAbelow(-0.5*dztilde) ...
        + dztilde^2*sRtilde(0);
    T = (Ktilde(dztilde) - zs*dztilde*(dKtilde(0) + vd)) ...
       /(Ktilde(-dztilde) + zs*dztilde*(dKtilde(0) + vd));
    Sb = 2*dztilde*Qtilde(xtilde) ...
        /(Ktilde(-dztilde) + zs*dztilde*(dKtilde(0) + vd));
    V = (Ktilde(ztilde(n-1)) + zs*dztilde*dKtilde(ztilde(n)))...
       /(Ktilde(ztilde(n)+dztilde) - zs*dztilde*dKtilde(ztilde(n)));
    Ctilde = zeros(n+1,m+1);
    Ctilde(:,1) = C0/Cs;
    for i=1:m
        A =  dxtilde*u0tilde(xtilde(i+1))*alphaTerm;
        E = -dxtilde*u0tilde( xtilde(i) )*alphaTerm;
        B = u0tilde(xtilde(i+1))*(4*dztilde^2*u0tilde( xtilde(i) ) + dxtilde*betaTerm);
        F = u0tilde( xtilde(i) )*(4*dztilde^2*u0tilde(xtilde(i+1)) - dxtilde*betaTerm);
        D =  dxtilde*u0tilde(xtilde(i+1))*gammaTerm;
        G = -dxtilde*u0tilde( xtilde(i) )*gammaTerm;
        
        Adiag = [nan,A(1) + T*D(1),A(2:end-1)];
        Ediag = [nan,E(1) + T*G(1),E(2:end-1)];
        Bdiag = B;
        Fdiag = F;
        Ddiag = [D(2:end-1),D(end) + V*A(end),nan];
        Gdiag = [G(2:end-1),G(end) + V*E(end),nan];
        sH = 2*dztilde^2*dxtilde*u0tilde(xtilde(i))*u0tilde(xtilde(i+1)) ...
            *(Stilde(xtilde(i),ztilde)+Stilde(xtilde(i+1),ztilde));
        
        sH(1) = sH(1) + D(1)*Sb(i)-G(1)*Sb(i+1);
    
        M = spdiags([Ddiag.', Bdiag.', Adiag.'],-1:1,n+1,n+1);
        N = spdiags([Gdiag.', Fdiag.', Ediag.'],-1:1,n+1,n+1);
        if n==20 && i==1
            full(M)
            full(N)
            full(sH.')
        end
        Ctilde(:,i+1) = N\(M*Ctilde(:,i) + sH.');
    end
end
clear variables
n = 20*2.^(0:7);
m = 5*2.^(0:7);

%%% Nondimensionalization Factors %%%
Cs = 1e6;
xs = 1e6;
zs = 1e6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Problem's Constants %%%
X = 10;
H = 200;
vd = 6e-3;
alpha_p = 0.61;
L = 1/0.0385;
zanchor = 100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

Kr = Kr_constant(alpha_p,L);

generate_conservative_form_functions(alpha_p,L,zanchor);

K = @(zvar) K_func(zvar,alpha_p,Kr);

schemeFunc = @myCDR;

myError = zeros(1,length(n));
C0 = 1;

Ctrue_func = @(xvar,zvar) C0 ...
    + xvar.*(H-zvar).^3;
syms xx zz
K0 = K_func(zanchor,alpha_p,Kr);
Cweight = sqrt(double(int(int((Ctrue_func(xs*xx,zs*zz+zanchor)/Cs).^2,zz,0,(H-zanchor)/zs),xx,0,X/xs)));
C_check_BC = K0*diff(Ctrue_func(xx,zz),zz);
C_check_BC = matlabFunction(C_check_BC);

if any(Ctrue_func(0,0:H) ~= C0) 
    Ctrue_func(0,0:H/10:H)
    error("Recheck your true solution to make sure it equals 0 at x=0")
elseif any(C_check_BC(0:X/10:X,H) ~= 0)
    error("Recheck your true solution to make sure the flux is 0 through z=H")
else
    Ctrue_func(0,0:H/10:H)
    C_check_BC(0:X/10:X,H)
    Q = @(xvar) vd*Ctrue_func(xvar,zanchor) - C_check_BC(xvar,zanchor);
end

mySource(xx,zz) = diff(Ctrue_func(xx,zz),xx) ...
    - (diff(sD(zz)*diff(Ctrue_func(xx,zz),zz) ...
    + sA(zz)*Ctrue_func(xx,zz),zz) + sR(zz)*Ctrue_func(xx,zz))./u0_func(xx);
S = matlabFunction(mySource);
Stilde = @(xvar,zvar) S(xs*xvar,zs*zvar+zanchor)*xs/Cs;

for i=1:length(n)
    Ctilde = schemeFunc(X,H,m(i),n(i),xs,zs,Cs,vd,alpha_p,Kr,C0,Q,zanchor,Stilde);
    dz = H/n(i);
    dx = X/m(i);
    [XX,ZZ] = meshgrid(0:X/m(i):X,0:H/n(i):H);
    Ctilde_true = Ctrue_func(xs*XX,zs*ZZ+zanchor);
    myError(i) = ...
        norm(Ctilde_true-Ctilde,"fro")*sqrt(dz/zs*dx/xs)/Cweight;
end

fig = figure;
loglog(1./n(1:end-0),myError(1:end-0),'-*b')
grid on
ylabel("Rel. Error","Interpreter","latex")
xlabel("$\Delta z$","Interpreter","latex")
myTitle = sprintf("$|| \\widetilde{C} - \\widetilde{C}_{true} ||_2\\:/\\:|| \\widetilde{C}_{true} ||_2$ with $\\alpha=%0.2f$",alpha_p);
tt = title(myTitle,"Interpreter","latex");

This produces the error plot:

Additionally, I use the following functions attached to run the setup.

[Torsten]

What is the PDE you are trying to solve ? What's the domain in which you want to solve it ? What are initial and boundary conditions ?

[David Gillcrist]

The PDE in nondimensionalized conservative form is:

The boundary conditions are

on the domain:

A lot more detail can be found in the linked Stack Exchange post on the problem itself. I didn't want to spend to much space here going through the same derivation when I could link it to another site I wrote it on already.

[Torsten]

Why not using "pdepe" to solve:

https://uk.mathworks.com/help/matlab/ref/pdepe.html

?

Of course, you must avoid u_0~ to become 0 for some value of x~ or z~.

[David Gillcrist]

I'll try to figure out how to use this function. How close 0 do you think $u_0$ can be before issues?

[Torsten]

For use of pdepe, you have to multiply your equation by u_0~ to get the derivative term "free" of a multiplying factor. And I'd define the flux f in "pdepe" as D*dC/dz, not as D*dC/dz + A*C. Instead, put the d/dz(A*C) together with S into the source term s in "pdepe". Otherwise, you will come into trouble when defining the boundary conditions.

How close 0 do you think $u_0$ can be before issues?

I usually set x~ and/or t~ to a small number to avoid division by 0 - it should be no problem to try which value(s) is/are appropriate and see whether it influences the result when the coding is finished.

[David Gillcrist]

@Torsten This ends up doing great! Thank you. Though I'm still concerned about why my scheme approach doesn't work.