Hi,
try:
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;
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
Best Answer