How to solve coupled odes with two time dependent variables with ode45?

I am having two coupled odes in which there are two dependent variables are present. I tried solving them with the below but I am getting error. Can someone help me in solving the below equations?

%%  equations which I need to solve are:
%% dc/dt = -((1-eps)*rhop/eps)*(dq/dt);
%% dq/dt = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n)) - q);

eps = 0.43;
rhop = 1228.5;
qm = 5.09;
R = 8.314;
T = 301;
Kl = 0.226;
n = 0.429;
K0 = 4.31e-9;           % in pascal-1
delh = -29380;          % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));

t = 0:1:100;
y0= zeros(2,1);
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);
plot(tsol,ysol(:,1))

function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
c = y(1);
q = y(2);
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
dy(2) = -((1-eps)*rhop/eps)*dy(1);

end

The error which I am getting when I am running this code is as follows:

Not enough input arguments.

Error in odcase>odfun (line 21)
dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);

Error in odcase>@(t,y)odfun(t,y) (line 14)
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in odcase (line 14)
[tsol,ysol] = ode45(@(t,y) odfun(t,y), t, y0);

Someone please let me know what are the mistakes in my code.

Answer 1

[Torsten]

function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
  q = y(1);
  c = y(2);
  dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
  dy(2) = -((1-eps)*rhop/eps)*dy(1);
end

instead of

function dy = odfun(t,y,Kl,qm,Keq,R,T,n,rhop,eps)
  c = y(1);
  q = y(2);
  dy(1) = Kl*((qm*Keq*c*R*T/(1+(Keq*c*R*T)^n)^(1/n))-q);
  dy(2) = -((1-eps)*rhop/eps)*dy(1);
end

Answer 2

[Sam Chak]

Hi @Ari Dillep Sai

Not sure what went wrong and why it is unstable. If you are absolutely sure that the absorption dynamics is stable (converging to a steady-state value), then the ODEs must be incorrect. Please check all parameters and the signs. Sometimes, a single change of sign can make a huge difference. For example, as a test, I simply added a minus sign '–' in front of K1 on the right-hand side of dydt(1), and the system becomes stable.

Please countercheck the ODEs against various textbooks and journal papers. If you only rely on a single reference and there is a misprint, then you know what happens...

function dydt = odefcn(t, y)

dydt = zeros(2,1);

ep = 0.43;

rhop = 1228.5;

qm = 5.09;

R = 8.314;

T = 301;

Kl = 0.226;

n = 0.429;

K0 = 4.31e-9; % in pascal-1

delh = -29380; % heat of adsorption in j/mol

Keq = K0*exp(-delh/(R*T));

c = y(1);

q = y(2);

dydt(1) = -Kl*( ( qm*Keq*c*R*T/(1 + (Keq*c*R*T)^n)^(1/n) ) - q );

dydt(2) = -((1 - ep)*rhop/ep)*dydt(1);

end

tspan  = [0 0.1];
init   = [1; 0];
[t, y] = ode45(@odefcn, tspan, init);
plot(t, y(:, 1), 'linewidth', 1.5)

[Sam Chak]

It appears that this question is related to another one here:

https://www.mathworks.com/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines?s_tid=srchtitle

and it looks like a boundary condition problem involving partial differential equations.

[Ari Dillep Sai]

Yeah, those are same equations. But I am not getting solution for those equations using method of lines. So first I thought of solving them by neglecting the spatial variations. So if I assume that there is no need of boundary conditions as they are varying only with z-axis.

[Ari Dillep Sai]

@Sam Chak Thank you for your reply!

I will once check the parameters and equations again.