% Script file to demonstrate poor long term performance of Matlab's ODE solvers
%
% Ian Mitchell for CS 403, 1/03/04

% ODE for a circular orbit.
f = inline('pi * [ -y(2); y(1) ]', 't', 'y');

% What solver?
odefun = @ode23;
%odefun = @ode15s;
%odefun = @ode23t;

% Get a phase plane plot during execution.
options = odeset('OutputFcn', @odeplot);


% Short time looks good.
[ ts, ys ] = feval(odefun, f, [ 0 4 ], [ 1; 0 ], options);

pause;

figure;
plot(ys(:,1), ys(:,2), 'b*-');
axis equal;
pause;


% Long time does not look good.
[ ts, ys ] = feval(odefun, f, [ 0 400 ], [ 1; 0 ], options);

pause;

figure;
plot(ys(:,1), ys(:,2), 'b*-');
axis equal;