matlab model integration

% simple model
clc;clear;close all
mu=0.5;k=100;
n0=1;
tend=30;
dt=0.1;
t=0:dt:tend;
nt=fix(tend/dt);
n=zeros(size(t));
n(1)=n0;
for i=1:length(t)-1
    dndt=mu*(1-n(i)/k)*n(i);
    n(i+1)=n(i)+dndt*dt;
end
plot(t,n)

　　

% model2
clc;clear;close all
mu1=0.5;s1=0.05;
mu2=0.4;s2=0.04;
n0=1;m0=1;
tend=30;dt=0.1;t=0:dt:tend;
nt=fix(tend/dt);
n=zeros(size(t));m=n;
n(1)=n0;m(1)=m0;
for i=1:length(t)-1
    dndt=(mu1-s1*(n(i)+m(i)))*n(i);
    dmdt=(mu2-s2*(n(i)+m(i)))*m(i);
    n(i+1)=n(i)+dndt*dt;
    m(i+1)=m(i)+dmdt*dt;
end
plot(t,n,t,m)

　　

 

 

% 2-3rd RK method:ode23
% 4-5th RK method:ode45
function dydt = modfun(t,y)
dydt =zeros(2,1);
mu1=0.5;s1=0.05;
mu2=0.4;s2=0.04;
dydt(1)=(mu1-s1*(y(1)+y(2)))*y(1);
dydt(2)=(mu2-s2*(y(1)+y(2)))*y(2);

　　

>> y0=[1,1];
>> [t,y]=ode23(@modfun,[0:30],y0);
>> plot(t,y)

　　

 

 

function dpdt = modfun(t,p)
dpdt =zeros(2,1);
fr=3.e13;pr=1.5;fx=60.e13;
v1=3.e16;v2=1.e18;
sed=0.01*fx*p(1);
dpdt(1)=(fr*pr-fx*p(1)+fx*p(2))/v1;
dpdt(2)=(fx*p(1)-fx*p(2)-sed)/v2;

　　

>> p0=[0,0]
>> [t,p]=ode23(@modfun,[0:1.e6],p0);
>> plot(t,p)

　

 

 

　

function dpdt = modfun(t,p)
dpdt =zeros(2,1);
fr=3.e13;pr=1.5;fx=60.e13;
v1=3.e16;v2=1.e18;tuo=1;
produc=p(1)*v1/tuo;
dpdt(1)=(fr*pr-fx*p(1)+fx*p(2)-produc)/v1;
dpdt(2)=(fx*p(1)-fx*p(2)+0.99*produc)/v2;

　　

>> [t,y]=ode23(@modfun,[0,1.e7],p0);
>> [t,y]=ode15s(@modfun,[0,1.e7],p);
>> plot(t,y)