function [ ] = prob_2_21
% 
% 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.
% 
%-----

y0 = [ 3; 0; 0 ];

options = odeset('OutputFcn',@odephas3,'Refine',1);
[ ts, ys ] = ode45( @f, [ 0 10 ], y0, options );

pause;

options = odeset('OutputFcn',@odephas3,'Refine',10);
[ ts, ys ] = ode45( @f, [ 0 10 ], y0, options );

return;



function [dydt] = f(t,y)
%
% y is 3x1
%
r = sqrt( y(1)^2 + y(2)^2 );
dydt = [ 
         -y(2) - y(1)*y(3)/r;
         +y(1) - y(2)*y(3)/r;
         +y(1)/r
       ];

return;