% solves delay differential equation
%
% dx/dt = alpha * x(t) * ( 1 - x(t-Tdelay) ) for t > 0, assuming
% x(t)=xs0 (constant) for -Tdelay <= t <0  
%
% method used is 2nd order RK
%
alpha=input('alpha= ')
Tdelay=1;
xs0=0.1;
%
% tmax is length of integration, dt is time step
%
tmax=20;
dt=0.01;
%
Ntmax=round(tmax/dt);
Ntdelay=round(1+Tdelay/dt);
xs=xs0*ones(1,Ntdelay);
xg=zeros(1,Ntmax+1);
xg(1)=xs0;
for it=1:Ntmax
    xn1=xs(Ntdelay)+alpha*xs(Ntdelay)*(1-xs(1))*dt;
    %xn2=0.5*xn1+0.5*(xs(Ntdel[ay)+alpha+xn1*(1-xs(1))*dt);
    xn2=xs(Ntdelay)+0.5*alpha*(xn1+xs(Ntdelay))*(1-0.5*(xs(1)+xs(2)))*dt;
    xs(1:Ntdelay-1)=xs(2:Ntdelay);
    xs(Ntdelay)=xn2;
    xg(it+1)=xn2;
end
hold off;times=[0:dt:tmax];plot(times,ones(1,Ntmax+1));axis([0 tmax 0 4]);hold on
%plot([0:dt:tmax],xg)
for it=1:Ntmax
    line([times(it) times(it+1)],[xg(it) xg(it+1)]);pause(0.0005);
end
hold off
