function sol = prob_4_14
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

x_b = -0.427; a = 0.16; epsilon = 0.02; u = 0.5; tau=1;

state = +1;
opts = ddeset('RelTol',1e-5,'Events',@events);
sol = dde23(@ddes,tau,[0.6],[0 120],opts,x_b,a,epsilon,u,tau,state);
while sol.x(end) < 120
    if sol.ie(end) == 1
        fprintf('New starting point at: %10.4f\n',sol.x(end));
        state = -state;
        sol = dde23(@ddes,tau,sol,[sol.x(end) 120],opts,x_b,a,epsilon,u,tau,state);
    else
        break;
    end
end
figure; plot(sol.x,sol.y(1,:),'-o'); grid on;

%=================================================================
function dydt = ddes(t,y,Z,x_b,a,epsilon,u,tau,state)
delta = Z(1)-x_b;
dydt = [ (1/tau)*( -y(1) + pi*( a + epsilon*state - u*(sin(delta))^2 ) ); ];

function [value,isterminal,direction] = events(t,y,Z,x_b,a,epsilon,u,tau,state)
value = [Z(1,1)-x_b];
isterminal = [1];
direction = [-state];