solving 13 interdependent odes simultaneously by ode45

qgin=12;
c_g_in=[1 2 3 4 5 6 ];
c=zeros(13,1);
c_l_in=[1 3 5 7 3 5 ];
vg=27;
kla=[13 34 5 6 7 7 ];
vl=65;
m=[23 3 5 2 1 4];
ql=123;
v_max=[1 2 3 5 7 1];
for i=1:6
    v(i)=-v_max(i).*c(i);
end
yxco=12;
R=8.314;
yxh2=23;
xp=45;
mu=-v(1)*yxco-v(2)*yxh2;
kd=15;
c0=[1 2 3 0 0 0 0 0 0 2 5 6 2];
c0g=c0(1:6);
psat=[1 23 45 54 67 86];
ngout=qgin*sum(c_g_in,'all')-kla(1)*((c(1)/m(1))-c(7))*vl-kla(2)*((c(2)/m(2))-c(8))*vl-kla(3)*((c(3)/m(3))-c(9))*vl+kla(4)*((c(10)/m(4))-c(4))*vl+kla(5)*((c(11)/m(5))-c(5))*vl+kla(6)*((c(12)/m(6))-c(6))*vl;
qgout=ngout*R*t/(sum((psat.*c0g),'all'));
[t,c]=ode45(@(t,c)odef(c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd),[0,600],c0)
function de = odef(c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd)
    de=zeros(1,13);   
de(1)=((qgin.*c_g_in(1)-qgout.*c(1))./vg)-kla(1).*(vl/vg).*((c(1)./m(1))-c(7));
de(2)=((qgin.*c_g_in(2)-qgout.*c(2))./vg)-kla(2).*(vl/vg).*((c(2)./m(2))-c(8));
de(3)=((qgin.*c_g_in(3)-qgout.*c(3))./vg)-kla(3).*(vl/vg).*((c(3)./m(3))-c(9));
de(4)=((qgin.*c_g_in(4)-qgout.*c(4))/vg)+kla(4).*(vl/vg).*((c(10)./m(4))-c(4));
de(5)=((qgin.*c_g_in(5)-qgout.*c(5))/vg)+kla(5).*(vl/vg).*((c(11)./m(5))-c(5));
de(6)=((qgin.*c_g_in(6)-qgout.*c(6))/vg)+kla(6).*(vl/vg).*((c(12)./m(6))-c(6));
de(7)=((ql/vl).*(c_l_in(1)-c(7))+kla(1).*((c(1)./m(1))-c(7))+v(1).*c(13));
de(8)=((ql/vl).*(c_l_in(2)-c(8))+kla(2).*((c(2)./m(2))-c(8))+v(2).*c(13));
de(9)=((ql/vl).*(c_l_in(3)-c(9))+kla(3).*((c(3)./m(3))-c(9))+v(3).*c(13));
de(10)=((ql/vl).*(c_l_in(4)-c(10))-kla(4).*((c(10)./m(4))-c4)+v(4).*c(13));
de(11)=((ql/vl).*(c_l_in(5)-c(11))-kla(5).*((c(11)./m(5))-c5)+v(5).*c(13));
de(12)=((ql/vl).*(c_l_in(6)-c(12))-kla(6).*((c(12)./m(6))-c6)+v(6).*c(13));
de(13)=((ql/vl).*(-c(13).*xp)+mu.*c(13)-kd.*c(13));

end
above is the code i wrote to solve the 13 odes simultaneously.i have done everything i could but i was unable to solve it. the errors kept popping out:

The following error occurred converting
from sym to double:
Unable to convert expression containing
symbolic variables into double array.
Apply 'subs' function first to
substitute values for variables.

Error in udhd>odef (line 29)
de(1)=((qgin.*c_g_in(1)-qgout.*c(1))./vg)-kla(1).*(vl/vg).*((c(1)./m(1))-c(7));

Error in
udhd>@(t,c)odef(c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd)
(line 26)
[t,c]=ode45(@(t,c)odef(c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd),[0,600],c0)

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 udhd (line 26)
[t,c]=ode45(@(t,c)odef(c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd),[0,600],c0)

the above code which i typed was in reference with this

    function dydt = odefcn(t,y,A,B)
    dydt = zeros(2,1);
    dydt(1) = y(2);
    dydt(2) = (A/B)*t.*y(1);
    A = 1;
    B = 2;
    tspan = [0 5];
    y0 = [0 0.01];
    [t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);

 the abive code i typed can be found on link below
<https://www.mathworks.com/help/matlab/ref/ode45.html solve non stiff odes with ode45>

Answer 1

[Alan Stevens]

Well, the following works. I wouldn't know if the results make sense!

qgin=12;
c_g_in=[1 2 3 4 5 6 ];
c=zeros(13,1);
c_l_in=[1 3 5 7 3 5 ];
vg=27;
kla=[13 34 5 6 7 7 ];
vl=65;
m=[23 3 5 2 1 4];
ql=123;
v_max=[1 2 3 5 7 1];
for i=1:6
    v(i)=-v_max(i).*c(i);
end
yxco=12;
R=8.314;
yxh2=23;
xp=45;
mu=-v(1)*yxco-v(2)*yxh2;
kd=15;
c0=[1 2 3 0 0 0 0 0 0 2 5 6 2];
c0g=c0(1:6);
psat=[1 23 45 54 67 86];
ngout=qgin*sum(c_g_in,'all')-kla(1)*((c(1)/m(1))-c(7))*vl-kla(2)*((c(2)/m(2))-c(8))*vl-kla(3)*((c(3)/m(3))-c(9))*vl+kla(4)*((c(10)/m(4))-c(4))*vl+kla(5)*((c(11)/m(5))-c(5))*vl+kla(6)*((c(12)/m(6))-c(6))*vl;
qgout=ngout*R/(sum((psat.*c0g),'all')); % t not yet defined so put it in odef
% ode45 must have t as first argument
% Need to pass v to odef
[t,c]=ode45(@(t,c)odef(t,c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd,v),[0,600],c0);

plot(t,c),grid
xlabel('t'),ylabel('c')

function de = odef(t,c,qgin,c_g_in,qgout,kla,vg,vl,m,c_l_in,ql,xp,mu,kd,v)

   qgout = qgout*t;  
    de=zeros(13,1);  % Needs to be a column vector   
de(1)=((qgin.*c_g_in(1)-qgout.*c(1))./vg)-kla(1).*(vl/vg).*((c(1)./m(1))-c(7));
de(2)=((qgin.*c_g_in(2)-qgout.*c(2))./vg)-kla(2).*(vl/vg).*((c(2)./m(2))-c(8));
de(3)=((qgin.*c_g_in(3)-qgout.*c(3))./vg)-kla(3).*(vl/vg).*((c(3)./m(3))-c(9));
de(4)=((qgin.*c_g_in(4)-qgout.*c(4))/vg)+kla(4).*(vl/vg).*((c(10)./m(4))-c(4));
de(5)=((qgin.*c_g_in(5)-qgout.*c(5))/vg)+kla(5).*(vl/vg).*((c(11)./m(5))-c(5));
de(6)=((qgin.*c_g_in(6)-qgout.*c(6))/vg)+kla(6).*(vl/vg).*((c(12)./m(6))-c(6));
de(7)=((ql/vl).*(c_l_in(1)-c(7))+kla(1).*((c(1)./m(1))-c(7))+v(1).*c(13));
de(8)=((ql/vl).*(c_l_in(2)-c(8))+kla(2).*((c(2)./m(2))-c(8))+v(2).*c(13));
de(9)=((ql/vl).*(c_l_in(3)-c(9))+kla(3).*((c(3)./m(3))-c(9))+v(3).*c(13));
de(10)=((ql/vl).*(c_l_in(4)-c(10))-kla(4).*((c(10)./m(4))-c(4))+v(4).*c(13));
de(11)=((ql/vl).*(c_l_in(5)-c(11))-kla(5).*((c(11)./m(5))-c(5))+v(5).*c(13));
de(12)=((ql/vl).*(c_l_in(6)-c(12))-kla(6).*((c(12)./m(6))-c(6))+v(6).*c(13));
de(13)=((ql/vl).*(-c(13).*xp)+mu.*c(13)-kd.*c(13));

end