ex21ch2.m

these codes are for solving OED with matlab

function ex21ch2
% torus problem
options = odeset('OutputFcn',@odephas3,'Refine',10);
ode45(@f,[0 10],[3; 0; 0],options);
%===================================================================
function dydt = f(t,y)
dydt = zeros(3,1);
r = sqrt(y(1)*y(1)+y(2)*y(2));
dydt(1) = -y(2) - (y(1) * y(3)) / r;
dydt(2) = y(1) - (y(2) * y(3)) / r;
dydt(3) = y(1) / r;