function prob_2_37
% 
% 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.
% 
%-----
  
n   = 3;
d   = 0.1;
yd  = 1 - (1/prod(1:3))*d^2 - (n/prod(1:5))*d^4;
ypd = - (2/prod(1:3))*d - (4*n/prod(1:5))*d^3;
erryd  = (5*n-8*n^2)*(d^6)/(3*(prod(1:7)));
errypd = 6*(5*n-8*n^2)*(d^5)/(3*(prod(1:7)));
fprintf('At d = %g, the error in y(d) is about %5.1e.\n',d,erryd);
fprintf('At d = %g, the error in y''(d) is about %5.1e.\n',d,errypd);

options = odeset('RelTol', 10^(-8), 'AbsTol', 10^(-10), 'Events',@events);

[x,y,xe,ye,ie] = ode45(@ode,[d 7],[ yd ypd ],options,n);
if( ~isempty(ie) ) 
  disp( 'Zero crossing at x=' ); 
  xe(1)
end

plot(x,y(:,1))

%========================================================
function dydx = ode(x,y,n)
  
dydx = [ y(2);
	 -(2/x)*y(2) - y(1)^n 
	 ];

%========================================================
function [value,isterminal,direction] = events(x,y,n)

value = y(1);
isterminal = 1; 
direction = 0;