% simulate Mackey-Glass equation
clear
%pick parameter values
delay = 6;
lambda = .05;
dt = .05;
k = fix(delay/dt);

%initialize x
x=ones(k,1);
X = [];


for j = 1:10000
xold = x(k);
xdel = x(1);
f= xdel.^3./(1+xdel.^3);

dxdt = lambda - xold.*f;
xnew = xold + dt*dxdt;
x = [x(2:k);xnew];

X(j) = xnew;
end
figure(1)
plot(X)


% plot neutral stabiity curve
gx= 2^(1/3);

y = [0:.01:gx];
f = y.^3./(1+y.^3);
fp = 3*y.^2./(1+y.^3).^2;
w= sqrt((y.*fp).^2-f.^2);


sigma = (pi + atan(-w./f))./w;
figure(2)
plot(y.*f,sigma)
axis([0 gx 0 20])




