function sol=minorsky
%
%   The dde is
%
%  y'=z;
%  z'=-y+epsilon*(-b1*z(t-tau)+c*z(t-tau).^3-z)
%

global tout
global tau
global b1
%tau = 3*pi;b1 = 3.8;
clf;
tau=7*pi; b1=3;



tout = 0;
options=ddeset('RelTol', 1e-5, 'AbsTol', 1e-5);


sol = dde23(@minde,[tau],@ddex1hist,[0, 400], options);
%figure;
clf;

plot(sol.x,sol.y(1,:));
hold on;
x=sol.x;
plot(x, 2*sqrt((1+b1*cos(tau))./3./cos(tau))+x*0, 'r');
% plot(sol.y(1,:), sol.y(2,:))
title(sprintf('b1=%g', b1));
xlabel('time t');
ylabel('solution y');

% --------------------------------------------------------------------------

function s = ddex1hist(t)
% Constant history function for DDEX1.
s = [1;0];

% --------------------------------------------------------------------------

function dydt = minde(t,Y,Z)
% Differential equations function for DDEX1.

global tau;
global b1;
epsilon = 0.1;
c = 1;


oo = Z(:,1);
ylag=oo(1);
zlag=oo(2);
y = Y(1);
z = Y(2);
dydt = [ z
         -y + epsilon*(-z-b1*zlag+c*zlag.^3)
      ];

