Slow integration of system ODE under specific conditions in adsorption problem

Question

I am trying to solve the system of ODEs below. I am applying the finite differences to z, so I get a system of ODE that I can solve with a solver like ode45.

The problem is that under the conditions that interest me integration is too slow (time step is very small) and I get no result. I would like to know what is the underlying problem and if there is any way to solve it. I have tried:

• using different solvers (tried all of them)
• using dimensionless variables
• using the pdepe function of Matlab to solve the PDE system directly without applying finite differences myself. Same problem occurs.

Full code (using the parameters that do not work)

function transport_model()
global parameter L Di epsb Q Cfeed K tfinal npz npt F ui h

parameter = [0.044511*432.75 432.75];    % isotherm parameters
L = 0.4;            % m, column length
epsb = 0.367;       % column bulk porosity
ui = 4e4;        % m/s, velocity
Cfeed = 200;     % feed concentration
K = 0.0782          % s-1, mass transfer coefficient
tfinal = 400;      % min, final time for calculation
npz = 50;           % number of discretization points in z
npt = 20;

F = (1-epsb)/epsb;
tspan = 0:tfinal/(npt-1):tfinal;
y0 = zeros(2*npz,1);
h = L/(npz-1);

sol = ode15s(@sedo, tspan, y0);
t = sol.x;
C1 = sol.y(1:npz,:);
q1 = sol.y((npz+1):end,:);

plot(t, C1(end,:))
xlim([0 tfinal])
ylim([0 Cfeed])



function DyDt = sedo(t,y)
global Cfeed K  npz F ui h

N = npz;
y(1) = Cfeed;
DyDt = zeros(2*N,1);

% Forward finite differences
DyDt(1) = -ui * 1/(2*h)*(-y(3)+4*y(2)-3*y(1)) - F * K * (y(1) - ilangmuir(y(N+1)));
DyDt(N+1) = K * (y(1) - ilangmuir(y(N+1)));

% Central finite differences
for i=2:N-1
    DyDt(i) = -ui * 1/(2*h)*(y(i+1)-y(i-1)) - F * K * (y(i) - ilangmuir(y(N+i)));
    DyDt(N+i) = K * (y(i) - ilangmuir(y(N+i)));
end

% Backward finite differences
DyDt(N) = -ui * 1/(2*h)*(3*y(N)-4*y(N-1)+y(N-2)) - F * K * (y(N) - ilangmuir(y(N+N)));
DyDt(2*N) = K * (y(N) - ilangmuir(y(N+N)));



function c = ilangmuir(q)
% langmuir solved for c
global parameter

a = parameter(1);
b = parameter(2);

c = q /(a - b * q);

However, if I try to model a different case (different conditions) I get the correct response. So I am assuming this is a numerical problem caused by the conditions I am trying to model and not an error in the model. For instance, if I use the conditions below I get the correct response.

parameter = [0.044511*432.75 432.75];    % isotherm parameters
L = 0.4;            % m, column length
epsb = 0.367;       % column bulk porosity
ui = 0.0056;        % m/s, velocity
Cfeed = 0.0025;     % feed concentration
K = 0.0782          % s-1, mass transfer coefficient
tfinal = 4000;      % min, final time for calculation
npz = 50;           % number of discretization points in z
npt = 20;