% this script considers the ODE y'(x) = x.y^2
% we aim to estimate y(1) given y(0)=1
% we do this using the Euler method with several different step sizes
% the Euler method requires the eulerFinal function 
% (a bit faster and more accurate than eulerSolve)

% this may seem complicated but a lot of the
% work is just plotting the graphs...

clc % clear screen
clear % clear all variables
format long

% parameters for the ODE
xMin = 0;
xMax = 1;
yInitial = 1;
Cval = 2/yInitial + xMin*xMin;

% choose some values of h to analyse
% for testing I suggest NN=10 
% ... for the data in the slides I used NN=16 
%     but it may take several hours to run
NN = 10;     % number of values
h(1) = 0.1;  % largest value of h
% each h is smaller than the last by a fixed ratio:
hRatio = 0.2;  
for i=2:NN
    h(i) = h(i-1)*hRatio;
end

% estimate y(xMax) for various h
for i=1:NN
    % we use the eulerFinal function, which estimates y(xMax) 
    yFinal(i) = eulerFinal(xMin,yInitial,xMax,h(i) );
end

% note yFinal(i) is our estimate of y(1), using step h(i)
% the following subtraction is performed for every element, to get the error
yErr = abs(2 - yFinal);
% plot the error as a function of h on a log-log plot
loglog(h,yErr,'bo-','Linewidth',1)
hold on
loglog(h,8*h,'k-')  % add comparison line
hold off
set(gca,'FontSize',18) 
set(gca, 'XTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
set(gca, 'YTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
legend('|E(1,h)| = |y(1)-Y(1,h)|','const * h','Location','NorthWest')
legend BOXOFF  % no box on legend
xlabel('h')
xlim( [ 1e-12 1] )
ylim( [ 1e-12 1] )

disp('** print error graph: eulerErr1.pdf');
% output to pdf
print -dpdf 'eulerErr1.pdf' -bestfit

% un-commment to explore further
%
% % solve the ODE again, using h values that are half as big
% for i=1:NN
%     h2 = 0.5*h(i);
%     yFinal2(i) = eulerFinal(xMin,yInitial,xMax,h2);
% end
% 
% % make a plot comparing step sizes h and h/2
% yDiff = abs( yFinal - yFinal2 );
% loglog(h,yDiff,'bo-','Linewidth',1)
% hold on
% loglog(h,4*h,'k-')
% hold off
% set(gca,'FontSize',18) 
% set(gca, 'XTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
% set(gca, 'YTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
% legend('| Y(1,h) - Y(1,h/2) | ','const * h','Location','NorthWest')
% legend BOXOFF  % no box on legend
% xlabel('h')
% xlim( [ 1e-12 1] )
% ylim( [ 1e-12 1] )
% 
% disp('** print stepsize comparison graph: eulerErr2.pdf');
% print -dpdf 'eulerErr2.pdf' -bestfit
% 
% % this is the "Richardson" estimate of y(x) at x=1
% yR = 2*yFinal2 - yFinal;
% % this is the error on the Richardson estimate
% yRerr = abs(2 - yR);
% % we plot only the first 10 points for the Richardson estimates
% % (after that the round-off errors get big)
% loglog(h(1:10),yRerr(1:10),'bo-','Linewidth',1)
% set(gca,'FontSize',18) 
% set(gca, 'XTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
% set(gca, 'YTick',[1e-12 1e-10 1e-8 1e-6 1e-4 1e-2 1])
% hold on
% loglog(h,yErr,'mo-','Linewidth',1)
% loglog(h,6*h,'k-')
% loglog(h,20*h.*h,'k--')
% legend('|Y_R(1,h)-y(1)|','|Y(1,h) - y(1)|','const * h',...
%     'const * h^2','Location','NorthWest')
% legend BOXOFF  % no box on legend
% hold off
% xlabel('h')
% xlim( [ 1e-12 1] )
% ylim( [ 1e-12 1] )
% 
% disp('** print Richardson error graph: eulerErr3.pdf');
% print -dpdf 'eulerErr3.pdf' -bestfit
% 
