%
% 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 = 50;
k = 0.05;
y_truth = @(t) ( y0*exp(-k*t) );

fprime = @(t,y) ( -k * y );

% Newton method: 
%
[tvals, yvals] = feuler( fprime, [0,10], y0, 1 ); v_1 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_1_newton = plot( tvals, relative_error, '-r' ); hold('on');

[tvals, yvals] = feuler( fprime, [0,10], y0, 0.1 ); v_2 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_2_newton = plot( tvals, relative_error, '-.r' );

[tvals, yvals] = feuler( fprime, [0,10], y0, 0.01 ); v_3 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_3_newton = plot( tvals, relative_error, ':r' );

[ v_1, v_2, v_3, y_truth(10) ]

% Trapezoidal method: 
%
[tvals, yvals] = eulertp( fprime, [0,10], y0, 1 ); v_1 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_1_trap = plot( tvals, relative_error, '-b' );

[tvals, yvals] = eulertp( fprime, [0,10], y0, 0.1 ); v_2 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_2_trap = plot( tvals, relative_error, '-.b' );

[ v_1, v_2, y_truth(10) ] 

% Runge-Kutta method: 
%
[tvals, yvals] = rkgen( fprime, [0,10], y0, 1, 1 ); v_1 = yvals(end); 
relative_error = ( yvals - y_truth(tvals) ) ./ y_truth(tvals);
ph_1_rk = plot( tvals, relative_error, '-g' ); 

[ v_1, y_truth(10) ] 

legend( [ph_1_newton,ph_2_newton,ph_3_newton,ph_1_trap,ph_2_trap,ph_1_rk], 'Euler h=1', 'Euler h=0.1', 'Euler h=0.01', 'Trapezoidal h=1', 'Trapezoidal h=0.1', 'RK: h=1', 'location', 'best' );
xlabel('time'); ylabel('y(t)');
%grid('on'); 

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