Ch4 in a Ch4/H2 breakthrough simulation occurs too quickly compared to experimental results.

I feel as though ive tried everything i can think of to fix it and so im here asking for suggestions as to what i should look into next. Any suggestions would be greatly appreciated. Many thanks

clear,clc
%% MAIN
% Physical Parameters
L   = 1.2;            % Column length m 
r   = 0.25;         % Bed radius m 

t = 0;
% tspan
t0  = 0;            % Initial Time
tf  = 2000;           % Final time
%dt  = 0.01;          % Time step
tspan   = [t0 tf];     % Time vector

dz  = 0.05;          % Mesh size
z   = 0:dz:L;       % Mesh generation
n   = numel(z);

%%%%%%% blended gas parameters %%%%%%%%%%%
%%% calc u input

TPF = [0.15 290 8.106e5 2e5 1e5];   % Feed: Superficial Velocity (m/s), tempurature (K) and pressure (Bar) [Vfeed Tfeed PH PI PL] 

yF  = 0.2;              % Mole fraction for methane 

%%%%%%% Preliminary calculations %%%%%%%%%%%

A = pi*r^2;         % Bed cross-sectional area m2
V = A*L;            % Bed volume m^3
mBed = A*L*426.7;   % Mass of adsorbent (in bed) 426.7 = bed density

%%%%%%%%% loop around here to work out p and n and mole fraction for each
%%%%%%%%% loop

% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
sol = adsorptionSolver(t,z,yF,TPF,tspan);

for i = 1:size(sol.y,2)
  [~,velocity(:,i),~,~] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF); 
end

for i = 1:size(sol.y,2)
  [~,~,vhalf(:,i),~] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);  
end

for i = 1:size(sol.y,2)
  [~,~,~,pspan(:,i),] = adsorptionModel(sol.x(i),sol.y(:,i),z(2)-z(1),numel(z),yF,TPF);  
end

% sol.x is time steps
% sol.y is 
Moley2 = 1-sol.y(1:n,1:end);

format long

purityh = 100*Moley2(n,end) / (sol.y(n,end) + Moley2(n,end));
fprintf('Purity of hydrogen is %f%% after %4.2fs\n', purityh, tf)


%% Balance Erros
% per meter^2 atm
epl = 0.433;
%MoleInCH4 = epl* 

dzhalf = dz/2;
zhalf   = (-dzhalf:dz:L+dzhalf)';
znodal = z';
zspan=cat(1,zhalf,znodal)';
zspan([1:2:end,2:2:end])=zspan; 

Cin = yF * TPF(3)/8.3145/TPF(2);
Cbed = sol.y(5*n,1:end).* sol.y(n,1:end)./8.3145./sol.y(4*n,1:end);
cratio = Cbed/Cin;

    
    figure(1)
    plot(sol.x,cratio)
    xlabel('C/C_0 of CH4 at bed end')
    ylabel('t')
    title('mole fraction ratio of methane')

%     figure(1)
%     
%     subplot(1,2,1)
%     mesh(sol.x,z,sol.y(1:n,1:end))
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('mole fraction y1')
%     title('mole fraction of methane')
%     
%     subplot(1,2,2)
%     mesh(sol.x,z,Moley2)
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('mole fraction y2')
%     title('mole fraction of hydrogen')
%     
%     figure(2)
%     
%     subplot(1,2,1)
%     mesh(sol.x,z,sol.y(n+1:2*n,1:end))
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('loading q1 mol/kg')
%     title('loading of methane')
%     
%     subplot(1,2,2)
%     mesh(sol.x,z,sol.y(2*n+1:3*n,1:end))
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('loading q2 mol/kg')
%     title('loading of hydrogen')
%     
%     figure(3)
%     
%     subplot(1,2,1)
%     mesh(sol.x,z,sol.y(3*n+1:4*n,1:end))
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('Temp')
%     title('temp of system')
%     subplot(1,2,2)
%     mesh(sol.x,z,sol.y(4*n+1:5*n,1:end))
%     xlabel('time')
%     ylabel('bed length')
%     zlabel('Pressure')
%     title('Pressure of system')
%     
%     figure(4)
%     subplot(1,2,1)
%     mesh(sol.x,z,velocity)
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('Velocity')
%     title('Intersistial velocity of system')
%     subplot(1,2,2)
%     mesh(sol.x,zspan,vhalf)
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('Velocity')
%     title('Intersistial velocity of system')
%     
%     
%     figure(5)
%     mesh(sol.x,z,pspan)
%     xlabel('t')
%     ylabel('bed length')
%     zlabel('dpdz')

   
%% Solver function
function out = adsorptionSolver(~,z,yFeed,TPFeed,tspan)



n = numel(z);      % Size of mesh grid
h = diff(z);
h = h(1);
% start up condidtions
Ta = 290;                   % Ambient Tempurature K
Pa = 8.106e5;                % Ambient Pressure Pa
y1a = 0;
y2a = 1-y1a;
Twa = 290;
%% LRC
Palrc = Pa./101325;
% CH4
k1c = 23.86e-3;	k2c = -0.0562e-3;	k3c = 2.81093e-3;	k4c = 1220; k5c = 1.628; k6c = -248.9;

% % H2
k1h = 7.34345e-3;	k2h = -0.013e-3;	k3h = 0.932e-3;	k4h = 506.306; k5h = 0.586972; k6h = 154.455;

qsi1 = k1c + k2c*Ta;
qsi2 = k1h + k2h*Ta;

B1 = k3c.*exp(k4c./Ta);
B2 = k3h.*exp(k4h./Ta);

n1 = (k5c+k6c./Ta);
n2 = (k5h+k6h./Ta); 

qstar1a = (qsi1.*B1.*(Palrc .*y1a).^n1) ./ (1+(((Palrc .*y1a).^n1.*B1)+((Palrc .*y2a).^n2.*B2)))*1e3; %mols/kg
qstar2a = (qsi2.*B2.*(Palrc .*y2a).^n2) ./ (1+(((Palrc .*y1a).^n1.*B1)+((Palrc .*y2a).^n2.*B2)))*1e3; %mols/kg

% Start up conditions Conditions
y1_0       = zeros(n,1) + y1a;     

q1_0        = zeros(n,1) + qstar1a;
q2_0        = zeros(n,1) + qstar2a;

T_0         = zeros(n,1) + Ta;

P_0         = zeros(n,1) + Pa;

Tw_0         = zeros(n,1) + Twa;
% y is an array size n*4 of c1 = 1:n, q1 = n+1:2*n, c2 = 2*n+1:3*n, q2 = 3*n+1:4*n
y0      = [y1_0; q1_0; q2_0; T_0; P_0; Tw_0];

out = ode15s(@(t,y) adsorptionModel(t,y,h,n,yFeed,TPFeed),tspan,y0);


end
%% Adsorption model
function [dydt, Velocity, Vhalf, Pspan] = adsorptionModel(t,y,h,n,yiFeed,TPFeed)
% set step no before run 
stepNo=1;


% Variables allocation
% y is an array size n*5 of y1 = 1:n, q1 = n+1:2*n,  
% q2 = 2*n+1:3*n, T = 3*n+1:4*n, P = 4*n+1:5*n
y1 = max(y(1:n),0);
y2 = 1 - y1;
y2 = max(y2,0);

q1 = max(y(n+1:2*n),0);
q2 = y(2*n+1:3*n);

T  = y(3*n+1:4*n);
P  = y(4*n+1:5*n);
Tw  = y(5*n+1:6*n);

% Z half points
yhalf = zeros(n+1,1);
Thalf = zeros(n+1,1);
Phalf = zeros(n+1,1);
uhalf = zeros(n+1,1);
deltPhalf = zeros(n+1,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Constant physical properties and basic simulation parameters
R = 8.3145;                 % ideal gas constant J/mol/K 
D   = 1.3e-5;               % Axial Dispersion coefficient m2/s
epl = 0.433;                 % Void fraction of bed
eplp = 0.61;               % Void fraction of particle
% eplt = epl + (1-epl)*eplp;  % Total void fraction
eplt = 0.758;  % Total void fraction 

% k   = [0.136;0.259];        % (lumped) Mass Transfer Coefficient 1/s
k   = [0.195;0.7];        % (lumped) Mass Transfer Coefficient 1/s
rhop = 850;               % Particle density kg/m3
% rhob = (1-epl)*rhop;
rhob = 532;
Cps = 1046.7;               % heat capacity of solid J/kg/K 
Kl = 0.09;                  % Thermal diffusivity J/m/s/K
% Kl = 1.2e-5;              % Thermal diffusivity J/m/s/K
deltaH = [23535; 12049.92];     % heat of adsorption J/mol
ufeed = TPFeed(1);
Tfeed = TPFeed(2);
PH = TPFeed(3); 
PI = TPFeed(4);
PL = TPFeed(5);


Tatm = 290;                   % ambient temp
hi = 38.4928;                   % heat transfer coefficient J/m2/s/K
ho = 14.2256;
% ignoring heat flux
% bed paramters
Rbi = 2.20447e-2;          % Bed inner radius m 
Rbo = 2.2073e-2;        % Bed outer raduis m
rhow = 7830;
Cpw = 502.416;
Aw = pi*(Rbo^2-Rbi^2);

% Ergun equation parameters
visc = 1.1e-5;     % gas viscosity kg/m/s
Rp = 1.15e-3;     % particle radius m

Acoeff = 1.75*(1-epl)/(2*Rp*epl^3);
Bcoeff = 150*(1-epl)^2/(4*Rp^2*epl^3);

% Constants for CH4
c1a = -0.703029; c1b = 108.4773; c1c = -42.52157; c1d = 5.862788; c1e = 0.678565;
% Constants for H2
c2a = 33.066178; c2b = -11.363417; c2c = 11.432816; c2d = -2.772874; c2e = -0.158558;

% Cpg = heat capacity (J/mol*K) 
Cpg1 = c1a + c1b*(T/1e3) + c1c*(T/1e3).^2 + c1d*(T/1e3).^3 + c1e./((T/1e3).^2);
Cpg2 = c2a + c2b*(T/1e3) + c2c*(T/1e3).^2 + c2d*(T/1e3).^3 + c2e./((T/1e3).^2);
Cpg = (y1.*Cpg1 + y2.*Cpg2);
% gas density
rhog = (P .* (y1*16.04e-3 + y2*2.02e-3)) ./ (R*T);  % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03

%%%%%%%% Boundary conditions: j = 0.5 and n + 0.5 %%%%%%%%%%%%%%%%%%%%%%%%    

   %% Adsorption (Ad)
if stepNo == 1
       uhalf(1) = ufeed - ufeed*exp(-0.5*t);
    %uhalf(1) = ufeed;
    
    Phalf(1) = P(1) + h/2*(Bcoeff*uhalf(1)*visc + Acoeff.*rhog(1)*uhalf(1)^2);
    
    %Phalf(1) = 15/8*P(1) -5/4*P(2) +3/8*P(3);
    %Phalf(1)  =  P(1) + (uhalf(1)*h*visc/2) / (4/150*(epl/(1-epl))^2*Rp^2)
    Phalf(n+1) = PH;

    yhalf(n+1) = y1(n);
    Thalf(n+1) = T(n);
    
    rhogOutlet = ((Phalf(n+1)) .* (yhalf(n+1)*16.04e-3 + (1-yhalf(n+1))*2.02e-3)) ./ (R*Thalf(n+1));
    rhogInlet = ((Phalf(1)) .* (yiFeed*16.04e-3 + (1-yiFeed)*2.02e-3)) ./ (R*Tfeed);
    Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
    Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
    CpgInlet = (yiFeed.*Cpg1in + (1-yiFeed).*Cpg2in);
    %
    
    deltPhalf(n+1) = (Phalf(n+1)-P(n));
    if deltPhalf(n+1) < 0
        uhalf(n+1) = (-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * 2/h*deltPhalf(n+1) .* (Acoeff.*rhogOutlet))) ./ (2*(Acoeff.*rhogOutlet));
    elseif deltPhalf(n+1) > 0
        deltPhalf(n+1) = -deltPhalf(n+1);
        uhalf(n+1) = -(-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * 2/h*deltPhalf(n+1) .* (Acoeff.*rhogOutlet))) ./ (2*(Acoeff.*rhogOutlet));
    else
        uhalf(n+1) = 0;
    end
    
    yhalf(1) = (y1(1)+yiFeed*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
    Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);

  %%  Purge (Pur)
elseif stepNo == 2
    
    uhalf(1) = ufeed - ufeed*exp(-0.5*t);
    
    Phalf(1) = P(1) + h/2*(Bcoeff*uhalf(1)*visc + Acoeff.*rhog(1)*uhalf(1)^2);
    Phalf(n+1) = PL;
    yhalf(n+1) = y1(n);
    Thalf(n+1) = T(n);

   
    %  yFeedfromAD = interp1(TimeoutAD,yAdOut,t);

    
    rhogOutlet = ((Phalf(n+1)) .* (yhalf(n+1)*16.04e-3 + (1-yhalf(n+1))*2.02e-3)) ./ (R*Thalf(n+1));
    rhogInlet = ((Phalf(1)) .* (yFeedfromAD*16.04e-3 + (1-yFeedfromAD)*2.02e-3)) ./ (R*Tfeed);
    Cpg1in = c1a + c1b*(Tfeed/1e3) + c1c*(Tfeed/1e3).^2 + c1d*(Tfeed/1e3).^3 + c1e./((Tfeed/1e3).^2);
    Cpg2in = c2a + c2b*(Tfeed/1e3) + c2c*(Tfeed/1e3).^2 + c2d*(Tfeed/1e3).^3 + c2e./((Tfeed/1e3).^2);
    CpgInlet = (yFeedfromAD.*Cpg1in + (1-yFeedfromAD).*Cpg2in);
    %
    
    deltPhalf(n+1) = (Phalf(n+1)-P(n));
    if deltPhalf(n+1) < 0
        uhalf(n+1) = (-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * 2/h*deltPhalf(n+1) .* (Acoeff.*rhogOutlet))) ./ (2*(Acoeff.*rhogOutlet));
    elseif deltPhalf(n+1) > 0
        deltPhalf(n+1) = -deltPhalf(n+1);
        uhalf(n+1) = -(-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * 2/h*deltPhalf(n+1) .* (Acoeff.*rhogOutlet))) ./ (2*(Acoeff.*rhogOutlet));
    else
        uhalf(n+1) = 0;
    end
    
    yhalf(1) = (y1(1)+yFeedfromAD*uhalf(1)*h/D/2) / (1+uhalf(1)*h/D/2);
    Thalf(1) = (T(1) + Tfeed*uhalf(1)*rhogInlet*CpgInlet*h/Kl/2)/(1 + uhalf(1)*rhogInlet*CpgInlet*h/Kl/2);   

end    
%% %%%%%% Wall values: j = 1.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) / (P(2)-P(1) + 1e-10)^4;
alpha1P = (1/3) / (2*(P(1)-Phalf(1)) + 1e-10)^4;
Phalf(2) = alpha0P/(alpha0P+alpha1P)*(0.5*(P(1)+P(2))) + alpha1P/(alpha0P+alpha1P)*(2*P(1)-Phalf(1));

alpha0y = (2/3) / (y1(2)-y1(1) + 1e-10)^4;
alpha1y = (1/3) / (2*(y1(1)-yhalf(1)) + 1e-10)^4;
yhalf(2) = alpha0y/(alpha0y+alpha1y)*(0.5*(y1(1)+y1(2))) + alpha1y/(alpha0y+alpha1y)*(2*y1(1)-yhalf(1));

alpha0T = (2/3) / (T(2)-T(1) + 1e-10)^4;
alpha1T = (1/3) / (2*(T(1)-Thalf(1)) + 1e-10)^4;
Thalf(2) = alpha0T/(alpha0T+alpha1T)*(0.5*(T(1)+T(2))) + alpha1T/(alpha0T+alpha1T)*(2*T(1)-Thalf(1));

%%%%%%%% Wall values: j = 2.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
alpha0P = (2/3) ./ (P(3:n)-P(2:n-1) + 1e-10).^4;
alpha1P = (1/3) ./ (P(2:n-1)-P(1:n-2) + 1e-10).^4;
Phalf(3:n) = alpha0P./(alpha0P+alpha1P).*(0.5.*(P(2:n-1)+P(3:n))) + alpha1P./(alpha0P+alpha1P).*(1.5.*P(2:n-1)- 0.5* P(1:n-2));

alpha0y = (2/3) ./ (y1(3:n)-y1(2:n-1) + 1e-10).^4;
alpha1y = (1/3) ./ (y1(2:n-1)-y1(1:n-2) + 1e-10).^4;
yhalf(3:n) = alpha0y./(alpha0y+alpha1y).*(0.5*(y1(2:n-1)+y1(3:n))) + alpha1y./(alpha0y+alpha1y).*(1.5.*y1(2:n-1)- 0.5* y1(1:n-2));

alpha0T = (2/3) ./ (T(3:n)-T(2:n-1) + 1e-10).^4;
alpha1T = (1/3) ./ (T(2:n-1)-T(1:n-2)+ 1e-10).^4;
Thalf(3:n) = alpha0T./(alpha0T+alpha1T).*(0.5*(T(2:n-1)+T(3:n))) + alpha1T./(alpha0T+alpha1T).*(1.5.*T(2:n-1)- 0.5* T(1:n-2));



%% %%%%%% Velocity wall values: j 1.5 to n - 0.5 %%%%%%%%%%%%%%%%%%%%%%%%
rhoghalf = (Phalf .* (yhalf*16.04e-3 + (1-yhalf)*2.02e-3)) ./ (R*Thalf); % ρ = PM/RT, M: H2 = 2.02 g/mol CH4 = 16.04 g/mol Air= 28.97 g/mol 50/50 H2 and CH4 = 9.03
deltPhalf(2:n) = (P(2:n) - P(1:n-1));
%deltPhalf(2:n) = (-Phalf(3:n+1)+8*P(2:n)-8*P(1:n-1)+Phalf(1:n-1))/6;

% Velocity calcs: uhalf(2:n) needs P(2:n) - P(1:n-1))
for i = 2:n
    if deltPhalf(i) < 0
        uhalf(i) = (-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * deltPhalf(i)/h .* (Acoeff.*rhoghalf(i)))) ./ (2*(Acoeff.*rhoghalf(i)));
    elseif deltPhalf(i) > 0
        deltPhalf(i) = -deltPhalf(i);
        uhalf(i) = -(-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * deltPhalf(i)/h .* (Acoeff.*rhoghalf(i)))) ./ (2*(Acoeff.*rhoghalf(i)));
    else
        uhalf(i) = 0;
    end
end

u = zeros(n,1);
deltP = (Phalf(2:n+1) - Phalf(1:n));

for i = 1:n
    if deltP(i) < 0
        u(i) = (-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * deltP(i)/h .* (Acoeff.*rhog(i)))) ./ (2*(Acoeff.*rhog(i)));
    elseif deltP(i) > 0
        deltP(i) = -deltP(i);
        u(i) = -(-(Bcoeff.*visc) + sqrt((Bcoeff.*visc)^2 - 4 * deltP(i)/h .* (Acoeff.*rhog(i)))) ./ (2*(Acoeff.*rhog(i)));
    else
        u(i) = 0;
    end
end


%% LRC
Plrc = P/101325;
% CH4
k1c = 23.86e-3;	k2c = -0.0562e-3;	k3c = 2.81093e-3;	k4c = 1220; k5c = 1.628; k6c = -248.9;

% % H2
k1h = 7.34345e-3;	k2h = -0.013e-3;	k3h = 0.932e-3;	k4h = 506.306; k5h = 0.586972; k6h = 154.455;

qsi1 = k1c + k2c*T;
qsi2 = k1h + k2h*T;

B1 = k3c.*exp(k4c./T);
B2 = k3h.*exp(k4h./T);

n1 = (k5c+k6c./T);
n2 = (k5h+k6h./T); 

qstar1 = (qsi1.*B1.*(Plrc.*y1).^n1) ./ (1+(((Plrc.*y1).^n1.*B1)+((Plrc.*y2).^n2.*B2)))*1e3; %mols/kg
qstar2 = (qsi2.*B2.*(Plrc.*y2).^n2) ./ (1+(((Plrc.*y1).^n1.*B1)+((Plrc.*y2).^n2.*B2)))*1e3; %mols/kg
LT1 = qstar1-q1;
LT2 = qstar2-q2;

%%%%%%%%%%%%%%% Eqautions for loading of component i %%%%%%%%%%%%%%
dq1dt = k(1)*LT1;
dq2dt = k(2)*LT2;

%%%%%%%%%%%%%% Balance equations

dPdt = zeros(n,1);
dy1dt = zeros(n,1);
dTdt = zeros(n,1);
dTwdt = zeros(n,1);

% Wall balance
dTwdt(1:n) = 1./(rhow*Cpw*Aw) .* ((2*pi*Rbi*hi*(T(1:n) - Tw(1:n)) - 2*pi*Rbo*ho*(T(1:n) - Tatm)));

% Energy balance
dTdt(1) = 1./(eplt.*rhog(1).*Cpg(1)+rhop*Cps) .* (Kl/h.*(((T(2)-T(1))/h)-2/h*(T(1)-Thalf(1))) - Cpg(1).*epl./h.*(rhoghalf(2).*Thalf(2).*uhalf(2)-rhoghalf(1).*Thalf(1).*uhalf(1)) + rhob.*(deltaH(1)*dq1dt(1)+deltaH(2)*dq2dt(1))  - 2*hi/Rbi*(T(1)-Tw(1)) );
dTdt(2:n-1) = 1./(eplt.*rhog(2:n-1).*Cpg(2:n-1)+rhop*Cps) .* (Kl/h.*(((T(3:n)-T(2:n-1))/h)-(T(2:n-1)-T(1:n-2))/h) - Cpg(2:n-1).*epl/h.*(rhoghalf(3:n).*Thalf(3:n).*uhalf(3:n)-rhoghalf(2:n-1).*Thalf(2:n-1).*uhalf(2:n-1)) + rhob.*(deltaH(1)*dq1dt(2:n-1)+deltaH(2)*dq2dt(2:n-1))  - 2*hi/Rbi*(T(2:n-1)-Tw(2:n-1)) );
dTdt(n) = 1./(eplt.*rhog(n).*Cpg(n)+rhop*Cps) .* (Kl/h.*(2/h*(Thalf(n+1)-T(n))-(T(n)-T(n-1))/h) - Cpg(n).*epl/h.*(rhoghalf(n+1).*Thalf(n+1).*uhalf(n+1)-rhoghalf(n).*Thalf(n).*uhalf(n)) + rhob.*(deltaH(1)*dq1dt(n)+deltaH(2)*dq2dt(n))  - 2*hi/Rbi*(T(n)-Tw(n)) );

% total mass balance
dPdt(1) = D.*T(1)/h.*( ((P(2)/T(2)-P(1)/T(1))/h)-2/h*(P(1)/T(1)-Phalf(1)/Thalf(1)))  -(T(1)./h).* (uhalf(2).*Phalf(2)./Thalf(2) - uhalf(1).*Phalf(1)./Thalf(1)) - ((rhop*R*T(1)*(1-epl))/epl).*(dq1dt(1)+dq2dt(1)) + P(1)./T(1).*dTdt(1);
dPdt(2:n-1) =D.*T(2:n-1)/h.*( (P(3:n)./T(3:n)-P(2:n-1)./T(2:n-1))/h - (P(2:n-1)./T(2:n-1)-P(1:n-2)./T(1:n-2))/h) -(T(2:n-1)./h).* (uhalf(3:n).*Phalf(3:n)./Thalf(3:n) - uhalf(2:n-1).*Phalf(2:n-1)./Thalf(2:n-1)) - ((rhop*R*T(2:n-1)*(1-epl))/epl).*(dq1dt(2:n-1)+dq2dt(2:n-1)) + P(2:n-1)./T(2:n-1).*dTdt(2:n-1);
dPdt(n) =D.*T(n)/h.*(2/h*(Phalf(n+1)./Thalf(n+1)-P(n)./T(n))-(P(n)./T(n)-P(n-1)./T(n-1))/h) -(T(n)./h).* (uhalf(n+1).*Phalf(n+1)./Thalf(n+1) - uhalf(n).*Phalf(n)./Thalf(n)) - ((rhop*R*T(n)*(1-epl))/epl).*(dq1dt(n)+dq2dt(n)) + P(n)./T(n).*dTdt(n);


% component mass balance
dy1dt(1) = (D.*T(1)./(h.*P(1))) * ( (Phalf(2)./Thalf(2)).*((y1(2)-y1(1))/h)-(Phalf(1)./Thalf(1)).*(2/h*(y1(1)-yhalf(1)))) - (T(1)./P(1)/h).*(uhalf(2).*Phalf(2).*yhalf(2)./Thalf(2) - uhalf(1).*Phalf(1).*yhalf(1)./Thalf(1)) - ((rhop*R*T(1)*(1-epl))/epl./P(1)).*(dq1dt(1)) - y1(1)./P(1).*dPdt(1) + y1(1)./T(1).*dTdt(1); 
dy1dt(2:n-1) = (D.*T(2:n-1)./(h.*P(2:n-1))) .* ( (Phalf(3:n)./Thalf(3:n)).*((y1(3:n)-y1(2:n-1))/h) - (Phalf(2:n-1)./Thalf(2:n-1)).*((y1(2:n-1)-y1(1:n-2))/h)) - (T(2:n-1)./P(2:n-1)/h).*(uhalf(3:n).*Phalf(3:n).*yhalf(3:n)./Thalf(3:n) - uhalf(2:n-1).*Phalf(2:n-1).*yhalf(2:n-1)./Thalf(2:n-1)) - ((rhop*R.*T(2:n-1).*(1-epl))./epl./P(2:n-1)).*(dq1dt(2:n-1)) - y1(2:n-1)./P(2:n-1).*dPdt(2:n-1) + y1(2:n-1)./T(2:n-1).*dTdt(2:n-1);
dy1dt(n) = (D.*T(n)./(h.*P(n))) * ( (Phalf(n+1)./Thalf(n+1)).*(2/h*(yhalf(n+1)-y1(n)))-(Phalf(n)./Thalf(n)).*((y1(n)-y1(n-1))/h)) - (T(n)./P(n)/h).*(uhalf(n+1).*Phalf(n+1).*yhalf(n+1)./Thalf(n+1) - uhalf(n).*Phalf(n).*yhalf(n)./Thalf(n)) - ((rhop*R*T(n)*(1-epl))/epl./P(n)).*(dq1dt(n)) - y1(n)./P(n).*dPdt(n) + y1(n)./T(n).*dTdt(n) ;% problem


%% caclutaions for graphs

Velocity = u;

Vhalf=cat(1,uhalf,u);
Vhalf([1:2:end,2:2:end])=Vhalf;

Pspan = (Phalf(2:n+1) - Phalf(1:n));

%%%%%%%%%%%%%%%% Concatenate vectors of time derivatives
dydt = [dy1dt;dq1dt;dq2dt;dTdt;dPdt;dTwdt];
end

[Torsten]

This seems to be a problem with your adsorption parameters (mass transfer coefficient, adsorption equilibrium, flow velocity ...). How should anyone out here be able to help in this case ?