4 non-linear differential equations with ODE45

Hi, i have a project about 4 non-linear differential equations where we have to show how temperature sex determination works in the reproductive systems of crocodiles. I have written a ode45 script/function which doesn't wnat to work/run. Is there someone who has had experience with ode45 and solving 4 different differential equation, if yes can you help me make my script/function work/run?

About the project:

Temperature-Dependent Sex Determination (TSD): Crocodilian Survivorship

• Literature: J.D. Murray. Mathematical biology: I. An introduction, volume 2. Springer, 2002. http://dl.icdst.org/pdfs/files/27f6eba850c27d335ff3f93778d8057f.pdf

-- In the chapter about TSD (page 144-155)

• A crucial difference between the crocodilia and most other species is that their sex is determined by the incubation temperature of the egg during gestation Because of TSD, crocodilia can recover from a catastrophic reduction in their population.

Model includes: females in region 1 (low temperature), females and males in region 2, males in region 3 (high temperature), leading to a system of 4 differential equations. • Question: how does environment (available place in regions) influence female-male ratio?

I have attached my script and function to this post.

Answer 1

[Stephan]

working scripts:

%%Vrdier 
% parametre:
% Carrying capacity
k1 = 500;
k2 = 600;
k3 = 700;
% De incubated males and females 
f1 = 500;
f2 = 400;
m2 = 400;
m3 = 500;
% birth rate
b = 1.5;
b0 = 1.5666;
% Death rate
d = 0.1;
% start values:
y0 = [10 20 10 20];
%time
time = [0 50];
%ode45 function
ty = @(t,x)crocofunction(t,x,b,k1,k2,k3,f1,f2,m2,m3,b0,d);
[t,y] = ode45(ty,time,y0);
% plot 
plot(t,y)

and:

function ty = crocofunction(t,x,b,k1,k2,k3,f1,f2,m2,m3,b0,d)
y0 = x(1);
y1 = x(2);
y2 = x(3);
y3 = x(4);
ty = zeros(4,1);
ty(1) = b*[k1./k1+f1]*f1-d*f1;
ty(2) = b0./2*[(f1^2./k1)+f2]*[k2./k2+f1+f2]-(d*f2);
ty(3) = b0./2*[(f1^2./k1)+f2]*[k2./k2+f1+f2]-(d*m2);
ty(4) = b0*[k3./k3+f1+f2]*[(f1^2./k1+f1)+f2]*[f1+f2./k2+f1+f2]-(d*m3);
end

but the result does not appear to be realistic. Are you using True world values? Then you should check your equations. From the view of pure Matlab now it works.