Solving set of time-dependent differential equations using ode45

I am writing a code to solve a set of two differential equations that are time-dependent, where μ is a constant:

I have written a code using a function odefcn as follows:

function dpddel = odefcn(p,del,mu)
dpddel = zeros(2,1);
dpddel(1) = (1/sqrt(mu))*sqrt(del(1-p)) - sqrt(mu)*p^(5/3);
dpdel(2) = sqrt(mu)*(1-p^(5/3)*del)/(1-p);
end

And then solving the set of equations using the ode45 solver:

% Defining constants
Co = 0.16;
H = 3.9;
S = 380;
cb = 0.6;
ct = 0.6;
ab = 4*1.095*0.25;
at = 4*2.125*2.8;
ab_star = sqrt(2)*cb*ab;
at_star = sqrt(2)*ct*at;
A_star = sqrt(1/((1/ab_star^2)+(1/at_star^2)));
B = 0.4325;
Td = (Co^(1/2)*S*H^(4/3))/(A_star*B^(1/3));
Tf = S/(Co*B^(1/3)*H^(2/3));
mu = Td/Tf;

% Solving equations
tspan = [0 60];
cond1 = p(0) == 1;
cond2 = del(0) == 1/Co;
conds = [cond1; cond2];
[p, del] = ode45(@(p,del) odefcn(p, del, mu), tspan, conds);
plot(p,del(:,1),'-o',p,del(:,2), '-.')

However the code will not solve the ode45 function, as it says there are not enough input arguments. I am not sure why this is the case. Is there a way of plotting the solutions to these two time-dependent differential equations?

Answer 1

[Stephan]

Hi,

try:

% Defining constants
Co = 0.16;
H = 3.9;
S = 380;
cb = 0.6;
ct = 0.6;
ab = 4*1.095*0.25;
at = 4*2.125*2.8;
ab_star = sqrt(2)*cb*ab;
at_star = sqrt(2)*ct*at;
A_star = sqrt(1/((1/ab_star^2)+(1/at_star^2)));
B = 0.4325;
Td = (Co^(1/2)*S*H^(4/3))/(A_star*B^(1/3));
Tf = S/(Co*B^(1/3)*H^(2/3));
mu = Td/Tf;

% Solving equations
tspan = [-10 60];
conds = [1; 1/Co];
[t, y] = ode45(@(t,x)odefcn(t,x,mu), tspan, conds);
plot(t,y(:,1),'-o',t,y(:,2), '-.')

function dpddel = odefcn(~,x, mu)
p = x(1);
del = x(2);
dpddel = zeros(2,1);
dpddel(1) = (1./sqrt(mu)).*sqrt(del*(1-p)) - sqrt(mu).*p.^(5/3);
dpdel(2) = sqrt(mu).*(1-p.^(5/3).*del)./(1-p);
end

Best regards

Stephan