Solve a pde equation with finite differences for Simulink

Hello , I want to transform this code that solves a pde equation with the ode solver into finite diferences, because I want to take the code as a matlab function block in simulink so it stands no ode solver(since it is an iterator take much time every time step so never ends simulation ) thats why i want to take it into finite differences .The equations are the following

The inital code is the following with ode solver:

L = 20    ;      % Longitud del lecho (m)
eps = 0.4;        % Porosidad
u = 0.2;       % Velocidad superficial del fluido (m/s)
k_f = 0.02;        % Constante de transferencia de masa (1/s)
c0 = 0;
Kf = 4;          % Constante de Freundlich
rhop = 1520;
n = 2;            % Exponente de Freundlich
         % Concentración inicial del fluido (kg/m³)
q0 = 4.320;        % Concentración inicial en el sólido (kg/m³)
      % Densidad del adsorbente (kg/m³)
tf = 10;         % Tiempo final de simulación (horas)
Nt = 100;

    t = linspace(0, tf*3600, Nt);
    Nz = 100;
    z = linspace(0, L,Nz);
    dz = z(2) - z(1);

    % Initial conditions
    ICA = max(ones(1, Nz) * c0, 1e-12); % Evitar valores negativos o cero
    ICB = ones(1, Nz) * q0;
    IC = [ICA ICB];
    options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'InitialStep', 1e-4, 'MaxStep', 100);

    [t, y] = ode15s(@fun_pde, t, IC, options, Nz, eps, n, Kf, k_f, u, rhop, dz);
    % Define value
    cc = y(:, 1:Nz);
    qq = y(:, Nz+1:end);
    % Recalculate new limit conditions
    cc(:, 1) = 0;
    cc(:, end) = cc(:, end-1);
    % Plotting
cp = cc(:, end) ./ c0;
qp = qq(:, :) ./ q0;
%q_promedio = mean(qq, 2);   % Promedio de q en el lecho para cada instante de tiempo
%conversion = 1 - (q_promedio / q0); % Conversión normalizada
figure;
subplot(2, 1, 1);
time = t / 3600; % Convertir a horas
plot(time, 1- qp, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Conversion');
title('Curva de conversión durante la desorción');
grid on;

subplot(2, 1, 2);
plot(t / 3600, (cc(:,:)), 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Soluciòn kg/m3');
title('Curva de carga de la solucion durante la desorciòn');
grid on;

    % PDE function
    function dydt = fun_pde(~, y, Nz, eps, n, Kf, k_f, u, rhop, dz)
        dcdt = zeros(Nz, 1);
        dqdt = zeros(Nz, 1);
        c = y(1:Nz);
        q = y(Nz+1:2*Nz);
        % Boundary conditions
        c(1) = max(c(1), 0); % Asegurar que c(1) sea no negativo
        c(end) = c(end-1); % Asegurar que c(1) sea no negativo
        % Interior nodes
        qstar = zeros(Nz, 1);
        dcdz = zeros(Nz, 1);
        for i = 2:Nz-1
            qstar(i) = Kf .* max(c(i), 1e-12).^(1/n); % Evitar problemas numéricos
            dqdt(i) = k_f .* (qstar(i) - q(i));
           % if i < Nz
                dcdz(i) = (c(i+1) - c(i-1)) / (2 * dz);
            %else
             %   dcdz(i) = (c(i) - c(i-1)) / dz;
            %end
            dcdt(i) = -u * dcdz(i) - rhop * ((1 - eps) / eps) .* dqdt(i);
        end
        dydt = [dcdt; dqdt];
    end
    
    
  

next is a try to solve with finite diferences but get someting different:

L = 20    ;      % Longitud del lecho (m)
eps = 0.4;        % Porosidad
u = 0.2;       % Velocidad superficial del fluido (m/s)
k_f = 0.02;        % Constante de transferencia de masa (1/s)
c0 = 0;
Kf = 4;          % Constante de Freundlich
rhop = 1520;
n = 2;            % Exponente de Freundlich
         % Concentración inicial del fluido (kg/m³)
q0 = 4.320;        % Concentración inicial en el sólido (kg/m³)
      % Densidad del adsorbente (kg/m³)
tf = 10;         % Tiempo final de simulación (horas)
Nz = 100;        % Número de nodos espaciales

% Discretización espacial y temporal
z = linspace(0, L, Nz);
t = linspace(0, tf*3600, Nt);
dz = z(2) - z(1);
dt = t(2) - t(1); % Paso temporal

% Condiciones iniciales
c = ones(Nt, Nz) * c0; % Concentración en el fluido
q = ones(Nt, Nz) * q0; % Concentración en el sólido

% Iteración en el tiempo (Diferencias Finitas Explícitas)
for ti = 1:Nt-1
    for zi = 2:Nz-1
        %  Isoterma de Freundlich
        qstar = Kf * max(c(ti, zi), 1e-12)^(1/n);

        % Transferencia de masa (Desorción)
        dqdt = k_f * (qstar - q(ti, zi));

        %  Gradiente espacial de concentración (Diferencias centradas)
        dcdz = (c(ti, zi+1) - c(ti, zi-1)) / (2 * dz);

        % Ecuación de balance de masa en el fluido
        dcdt = -u * dcdz - rhop * ((1 - eps) / eps) * dqdt;

        % Actualizar valores asegurando que sean positivos
        c(ti+1, zi) = max(c(ti, zi) + dcdt * dt, 0);
        q(ti+1, zi) = max(q(ti, zi) + dqdt * dt, 0);
    end
end

% Condiciones de frontera
c(:, 1) = c0;        % Entrada con concentración baja
c(:, Nz) = c(:, Nz-1); % Gradiente nulo en la salida
% Cálculo de la conversión normalizada
qp = q(:, :) ./ q0;


% Graficar resultados
figure;
subplot(2, 1, 1);
plot(t / 3600, 1-qp, 'b', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Conversion');
title('Curva de conversión durante la desorción');
grid on;

subplot(2, 1, 2);
c_salida = c(:, :); % Concentración en la salida del lecho
plot(t / 3600, c_salida, 'r', 'LineWidth', 1.5);
xlabel('Tiempo (horas)');
ylabel('Soluciòn kg/m3');
title('Curva de carga de la solucion durante la desorciòn');
grid on;

I dont know why is wrong

Thanks in advance