MATLAB: System of differential equations with many input arguments.

differential equationsMATLABMATLAB and Simulink Student Suiteode45odefcn

Hello,
I am trying to solve a system of five differential equations and get some plots back. The equations I have to solve are the following:
Where all correspond to derivatives of specific functions which in order to make the code for matlab I substituted them in the above functions. The code I am using for that is the following:
% l stands for \tau, u for \Theta, v for a, o for b (from utility) and j for ξ
syms S(t) P(t) c(t) h(t) z(t) A b d e g k l m n o r t u v w Y
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - v*z - c - b*S*l;
eqn2 = diff(P,t) == u*z - d*P + j*c + (1-b)*S;
eqn3 = diff(h,t) == m*((b*S)-h);
eqn4 = diff(c,t) == (((b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l -r)) - ((A*g*(b*S)^(1-g)*z^(g-1) - v) * (1 - b + j*(d+r)) * (1/u))) * (u/(u - j*(v-g*A*(b*S)^(1-g)*z^(g-1))))) + ((1/(n*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) * ((exp(r*t)*w*m*b + b*(c*P^(-e)*h^(-o))^(1-n)*S^(-n)*b^(1-n)) + (j*e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)) - ((1-n)*e*c^(-n)*P^(-e*(1-n)-1)*(b*S*h^(-o))^(1-n)*(u*z+(1-b)*S+j*c-d*P)) + (b*(1-n)*c^(-n)*(P^(-e)*h^(-o))^(1-n)*b^(1-n)*S^(-n)*(A*(b*S)^(1-g)*z^g - c - v*z - l*b*S))));
eqn5 = diff(z,t) == (((A*g*(b*S)^(1-g)*z^(g-1)-v) / (A*g*(g-1)*(b*S)^(1-g)*z^(g-2))) * ( d + ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*u) / (c^(-n)*(P*b*S*h^(-o))^(1-n))) - ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*j) / (c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l + ((w*m*(u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (exp(-r*t)*u*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n)))) + ((b^(1-n)*S^(-n)*(c*P^(-e)*h^(-o))^(1-n) * (u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (u*e^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + ((1-b)*(v-A*g*(b*S)^(1-g)*z^(g-1))*(1/u)))) - ((b*z/S) * (A*(b*S)^(1-g)*z^g - v*z - c - b*S*l));
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4,eqn5);
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w});
A=3;
b=0.6;
d=0.3;
e=0.4;
g=0.3;
k=4;
l=2;
m=0.7;
n=2;
o=0.6;
r=0.5;
u=10;
v=2;
w=0.7;
tspan = [0 50];
y0 = [100 60 40 30 50];
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w) , tspan, y0);
P = Y(:,1);
S = Y(:,2);
h = Y(:,3);
c = Y(:,4);
z = Y(:,5);
x = b*S;
q = A*((b*S).^(1-g)).*(z.^g);
figure(1)
plot(c,P)
title(sprintf('c-P'))
xlabel('Consumption (c)')
ylabel('Pollution (P)')
figure(2)
plot(c,x)
title(sprintf('C-x'))
xlabel('Consumption (c)')
ylabel('Recycling (x)')
figure(3)
plot(q,P)
title(sprintf('q-P'))
xlabel('Output (q)')
ylabel('Pollution (P)')
figure(4)
plot(q,x)
title(sprintf('q-x'))
xlabel('Output (q)')
ylabel('Recycling (x)')
I've used that code before for different functions and it was working perfectly but with those differential equations it doesn't work and it give me back the following error:
Error using
symengine>@(t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w)[-d.*Y(1)+u.*Y(5)+Y(4).*1i-Y(2).*(b-1.0);-v.*Y(5)-Y(4)+A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g-b.*l.*Y(2);m.*(b.*Y(2)-Y(3));-(u.*(b.*(l+r+A.*b.^(-g+1.0).*Y(2).^(-g).*Y(5).^g.*(g-1.0))-((v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(-b+d.*1i+r.*1i+1.0))./u))./(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i)+(Y(4).^n.*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0).*(e.*Y(1).^(e.*(n-1.0)-1.0).*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0).*1i+b.*m.*w.*exp(r.*t)+b.*b.^(-n+1.0).*Y(2).^(-n).*(Y(1).^(-e).*Y(3).^(-o).*Y(4)).^(-n+1.0)-e.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^(-n).*(n-1.0).*(b.*Y(2).*Y(3).^(-o)).^(-n+1.0).*(d.*Y(1)-u.*Y(5)-Y(4).*1i+Y(2).*(b-1.0))+b.*b.^(-n+1.0).*(Y(1).^(-e).*Y(3).^(-o)).^(-n+1.0).*Y(2).^(-n).*Y(4).^(-n).*(n-1.0).*(v.*Y(5)+Y(4)-A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g+b.*l.*Y(2))))./n;(b.*Y(5).*(v.*Y(5)+Y(4)-A.*(b.*Y(2)).^(-g+1.0).*Y(5).^g+b.*l.*Y(2)))./Y(2)-((b.*Y(2)).^(g-1.0).*Y(5).^(-g+2.0).*(v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(d-b.*(l+A.*b.^(-g+1.0).*Y(2).^(-g).*Y(5).^g.*(g-1.0)-(m.*w.*exp(r.*t).*Y(4).^n.*(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0))./u)-((v-A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0)).*(b-1.0))./u-e.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^n.*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0).*1i+e.*u.*Y(1).^(e.*(n-1.0)-1.0).*Y(4).^n.*(b.*Y(1).*Y(2).*Y(3).^(-o)).^(n-1.0).*(b.*Y(2).*Y(3).^(-o).*Y(4)).^(-n+1.0)+(b.^(-n+1.0).*e.^n.*Y(2).^(-n).*(u-v.*1i+A.*g.*(b.*Y(2)).^(-g+1.0).*Y(5).^(g-1.0).*1i).*(Y(1).^(-e).*Y(3).^(-o).*Y(4)).^(-n+1.0).*(b.*Y(1).^(-e).*Y(2).*Y(3).^(-o)).^(n-1.0))./u))./(A.*g.*(g-1.0))]
Too many input arguments.
Error in @(t,Y)odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w)
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);
So now I do not know how to move on. I understand that the system of differential equations I am using has many arguments as the error message says as well, but most of them are parameters which I define giving numerical values for them. Shouldn't Matlab be able to solve that?
Any help would be very much appreciated (any idea, thought, solution or anything else).
Thank you so much in advance.

Best Answer

Try:
% l stands for \tau, u for \Theta, v for a, o for b (from utility) and j for ξ
syms S(t) P(t) c(t) h(t) z(t) A b d e g k l m n o r t u v w Y
eqn1 = diff(S,t) == A*(b*S)^(1-g)*(z^g) - v*z - c - b*S*l;
eqn2 = diff(P,t) == u*z - d*P + j*c + (1-b)*S;
eqn3 = diff(h,t) == m*((b*S)-h);
eqn4 = diff(c,t) == (((b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l -r)) - ((A*g*(b*S)^(1-g)*z^(g-1) - v) * (1 - b + j*(d+r)) * (1/u))) * (u/(u - j*(v-g*A*(b*S)^(1-g)*z^(g-1))))) + ((1/(n*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) * ((exp(r*t)*w*m*b + b*(c*P^(-e)*h^(-o))^(1-n)*S^(-n)*b^(1-n)) + (j*e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)) - ((1-n)*e*c^(-n)*P^(-e*(1-n)-1)*(b*S*h^(-o))^(1-n)*(u*z+(1-b)*S+j*c-d*P)) + (b*(1-n)*c^(-n)*(P^(-e)*h^(-o))^(1-n)*b^(1-n)*S^(-n)*(A*(b*S)^(1-g)*z^g - c - v*z - l*b*S))));
eqn5 = diff(z,t) == (((A*g*(b*S)^(1-g)*z^(g-1)-v) / (A*g*(g-1)*(b*S)^(1-g)*z^(g-2))) * ( d + ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*u) / (c^(-n)*(P*b*S*h^(-o))^(1-n))) - ((e*P^(-e*(1-n)-1)*(c*b*S*h^(-o))^(1-n)*j) / (c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + b*(A*(1-g)*b^(1-g)*S^(-g)*z^g - l + ((w*m*(u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (exp(-r*t)*u*c^(-n)*(P^(-e)*b*S*h^(-o))^(1-n)))) + ((b^(1-n)*S^(-n)*(c*P^(-e)*h^(-o))^(1-n) * (u-j*(v-A*g*(b*S)^(1-g)*z^(g-1)))) / (u*e^(-n)*(P^(-e)*b*S*h^(-o))^(1-n))) + ((1-b)*(v-A*g*(b*S)^(1-g)*z^(g-1))*(1/u)))) - ((b*z/S) * (A*(b*S)^(1-g)*z^g - v*z - c - b*S*l));
[VF,Subs] = odeToVectorField(eqn1,eqn2,eqn3,eqn4,eqn5);
odefcn = matlabFunction(VF, 'Vars',{t,Y,A,b,d,e,g,k,l,m,n,o,r,u,v,w});
A=3;
b=0.6;
d=0.3;
e=0.4;
g=0.3;
k=4;
l=2;
m=0.7;
n=2;
o=0.6;
r=0.5;
u=10;
v=2;
w=0.7;
tspan = [0 50];
y0 = [100 60 40 30 50];
[t,Y] = ode45(@(t,Y) odefcn(t,Y,A,b,d,e,g,k,l,m,n,o,r,u,v,w) , tspan, y0);
P = Y(:,1);
S = Y(:,2);
h = Y(:,3);
c = Y(:,4);
z = Y(:,5);
x = b*S;
q = A*((b*S).^(1-g)).*(z.^g);
figure(1)
plot(c,P)
title(sprintf('c-P'))
xlabel('Consumption (c)')
ylabel('Pollution (P)')
figure(2)
plot(c,x)
title(sprintf('C-x'))
xlabel('Consumption (c)')
ylabel('Recycling (x)')
figure(3)
plot(q,P)
title(sprintf('q-P'))
xlabel('Output (q)')
ylabel('Pollution (P)')
figure(4)
plot(q,x)
title(sprintf('q-x'))
xlabel('Output (q)')
ylabel('Recycling (x)')
3 reasons:
[t,Y] = ode45(@(t,Y) odefcn(t,A,Y,b,d,e,g,k,l,m,n,o,r,t,u,v,w) , tspan, y0);
% ^ ^ ^
% | | |
% -----> Change order ----> t already at first position
odefcn = matlabFunction(VF, 'Vars',{t,A,Y,b,d,e,g,k,l,m,n,o,r,u,v,w});
% ^ ^
% | |
% SAME HERE -----> Change order