function [ ] = prob_2_26
%
% Note: I couldn't get this to work (because matlab would crash when
% I had more than about 4 elements in xarray.  So I created a for loop 
% for each.
% 
% Written by:
% -- 
% John L. Weatherwax                2005-04-03
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

global xindx xarray;
xarray = 0.1:0.1:0.9;

close all; 

options = odeset('Events',@events,'RelTol',1d-5,'AbsTol',1d-10);


y0 = [ 0; 0 ]; 
ter=[];
for xi=1:length(xarray)
  xindx=xi; 
  [ts,ys,te,ye,ie] = ode45( @f0, [0, 2], y0, options );
  ter = [ ter, te ];
end
figure; plot( ts, ys(:,1), '-or' ); title( 'dydx=f0 solution' ); 

figure;
plot( xarray, ter, '-or' ); hold on;
plot( xarray, 2*sqrt(xarray), '-+g' ); grid on; title( 'event times and exact' ); 
figure; plot( xarray, ter-2*sqrt(xarray), '-ob' ); title( 'approx-exact' ); 
drawnow;


y0 = [ 0; 1 ]; 
ter=[];
for xi=1:length(xarray)
  xindx=xi; 
  [ts,ys,te,ye,ie] = ode45( @f1, [0, 2], y0, options );
  ter = [ ter, te ];
end
figure; plot( ts, ys(:,1), '-or' ); title( 'dydx=f1 solution' ); 

figure;
plot( xarray, ter, '-or' ); hold on;
plot( xarray, 2-2*sqrt(1-xarray), '-+g' ); grid on; title( 'event times and exact' ); 
figure; plot( xarray, ter-2*(1-sqrt(1-xarray)), '-ob' ); title( 'approx-exact' ); 
drawnow;


y0 = [ 0; 1 ]; 
ter=[];
for xi=1:length(xarray)
  xindx=xi; 
  [ts,ys,te,ye,ie] = ode45( @f2, [0, 2], y0, options );
  ter = [ ter, te ];
end
figure; plot( ts, ys(:,1), '-or' ); title( 'dydx=f2 solution' ); 
  
figure;
plot( xarray, ter, '-or' ); hold on;
plot( xarray, asin(xarray), '-+g' ); grid on; title( 'event times and exact' ); 
figure; plot( xarray, ter-(asin(xarray)), '-ob' ); title( 'approx-exact' ); 
drawnow;  

return;



function [dydt] = f0(t,y)
%
% y is 2x1
%
dydt = [ 
         y(2);
         0.5
       ];

return;

function [dydt] = f1(t,y)
%
% y is 2x1
%
dydt = [ 
         y(2);
         -0.5
       ];

return;

function [dydt] = f2(t,y)
%
% y is 2x1
%
dydt = [ 
         y(2);
         -y(1)
       ];

return;



function [value,isterminal,direction] = events(t,y)

global xindx xarray;

% Old code (didn't work when xarray was an entire vector?):
% $$$ value = [ xarray - y(1) ].';
% $$$ isterminal = zeros(length(xarray),1);
% $$$ direction = zeros(length(xarray),1);

% New code:
value = [ xarray(xindx) - y(1) ].';
isterminal = 0; 
direction = 0;

return;