function do_all=delay_Schopf_Suarez_1989
% A matlab version of the model of Schopf and Suarez 1989.
% Eli Tziperman, for EPS231

close all;
clear all;

% set parameters for the delay equation
DT=0.1;  % time is in months
alpha=0.54;  % 0.75: standard value; 0.54: growing oscillations; 0.51: decaying
delta=round(8/DT);

% set initial conditions:
i_start=round(delta+1);
i_end=round(1000/DT);
% start at fixed point plus a small perturbation
h(1:i_start-1)=sqrt(1-alpha)+0.001; 
h(i_start:i_end)=0.0;

% integrate the delay equation in time:
for ii=i_start:i_end
  h(ii)=h(ii-1)+DT*(h(ii-1)-alpha*h(ii-delta)-h(ii-1)^3);
end

% time axis in years for plotting:
time=[DT/12:DT/12:DT*i_end/12];

figure
% plot time series:
subplot(2,1,1)
is=i_start;ie=i_end;
plot(time(is:ie),h(is:ie))
h1=xlabel('time (years)');
h2=ylabel('h');
h3=title('time series');
set([gca,h1,h2,h3],'fontsize',20)

% plot reconstructed delay coordinates phase space:
subplot(2,1,2)
is=round(i_start+2*delta);
ie=i_end;
plot(h(is-2*delta:ie-2*delta),h(is:ie),'-')
h1=xlabel('h(t)');
h2=ylabel('h(t-year)');
h3=title('reconstructed phase space');
set([gca,h1,h2,h3],'fontsize',20)
