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 advanceHello , 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