clear all; close all;
% simulate Mackey_Glass equations
% udot(t)=-gamma*u(t)+beta*u(t-tau)/(1+u(t-tau)^n)

% parameters
gamma=1;
beta=2;
n=20;
tau=2;

% Mackey-Glass function
mgfun=@(t,y,Z) -gamma*y+beta*Z/(1+Z^n);

% specify a constant history not equal to either steady state
history=2*(beta/gamma-1)^(1/n);

% sufficiently long interval to evole through transient
tspan=[0 1000];

% use enough plot points to see smooth solution
plotpoints=20*(tspan(2)-tspan(1));

% using enough plotpoints produces a nearly smooth solution
% we also need to reduce the tolerances from their default values
% options = ddeset('AbsTol',1e-6,'RelTol',1e-3);
options = ddeset('AbsTol',1e-8,'RelTol',1e-6);

% use matlab dde23 solver to solve ivp
sol = dde23(mgfun,tau,history,tspan,options);

tint=linspace(tspan(1),tspan(2),plotpoints);
uval = deval(sol,tint);
figure(1)
plot(tint,uval)

figure(2)
tint=linspace(0.2*(tspan(1)+tau)+0.8*tspan(2),tspan(2),plotpoints);
uval = deval(sol,tint);
tdelint=tint-tau;
udelval = deval(sol,tdelint);
plot(uval,udelval);

fil=['PPdde23-',num2str(n),'.eps'];
figfil=['PPdde23-',num2str(n),'.fig'];
print('-depsc',fil)
saveas(gcf,figfil,'fig')