%
% Written by:
% --
% John L. Weatherwax                2009-04-21
%
% email: wax@alum.mit.edu
%
% Please send comments and especially bug reports to the
% above email address.
%
%-----

close all; clc; clear all; 

addpath('../../BookCode');

y0 = [1/2; -1/2];

y_truth = @(t) ( 1 ./ ( 1 + t.^2 ) );

fprime = @(t,y) ( [ y(2); y(2)/t + 8 * t^2 * y(1)^3 ] );

[tvals, yvals] = ode23( fprime, [1,4], y0 ); y_approx = yvals(:,1); v_1 = yvals(end,1); 
relative_error = ( y_approx - y_truth(tvals) ) ./ y_truth(tvals);
ph_ode23 = plot( tvals, relative_error, '-r' ); hold('on');

[tvals, yvals] = ode45( fprime, [1,4], y0 ); y_approx = yvals(:,1); v_2 = yvals(end,1); 
relative_error = ( y_approx - y_truth(tvals) ) ./ y_truth(tvals);
ph_ode45 = plot( tvals, relative_error, '-b' ); 

[ v_1, v_2, y_truth(4) ] 

legend( [ph_ode23,ph_ode45], 'ode23', 'ode45', 'location', 'southwest' );
xlabel('time'); ylabel('y(t)');
grid('on'); 

%saveas( gcf, '../../WriteUp/Graphics/Chapter5/c_5_p_4_plot.eps', 'epsc' );