% script: exer016Test.m
% description: solves exercise0 of NumAnal Seminar 2016: 
% dy/dt = -t*y*y,  0<=t<=3
% y(0) = 1
% solution y(t) = 2/(t^2 + 2)

 % programmer: t. shores
% last rev: 1/30/16

% Octave and Matlab compatible file
disp('Exercise 0, NumAnal Seminar 2016:')
disp('dy/dt = -t*y*y,  0<=t<=3')
disp('y(0) = 1')
% initialization code
t0 = 0;
tLast = 3;
y0 = exer016Soln(t0);
fun = 'exer016RHS';
numSteps = 8*[1,2,4,8];
hSteps = (tLast - t0)./numSteps;
lenSteps = length(hSteps);
errors = zeros(lenSteps, 4);
% store stepsizes and errors in this matrix
errors(:,1) = hSteps;
% main  -- we use the infinity norm for error
ySoln = exer016Soln(tLast);
for jj = 1:lenSteps
  yApprox = Euler(fun,t0,y0,hSteps(jj),numSteps(jj));
  errors(jj,2) = norm(yApprox - ySoln,inf);
  yApprox = Trap(fun,t0,y0,hSteps(jj),numSteps(jj));
  errors(jj,3) = norm(yApprox - ySoln,inf);
  yApprox = ClassERK(fun,t0,y0,hSteps(jj),numSteps(jj));
  errors(jj,4) = norm(yApprox - ySoln,inf);
end
disp('Table with columns:')
disp(' [Time stepsize, Euler error, trapezoid error, ERK4 error]');
disp(errors)
disp('Table with columns of rates of decrease of error');
disp(' [Stepsize,   Euler,   trapezoid,   ERK4]');
disp(errors(1:lenSteps-1, 1:4)./errors(2:lenSteps, 1:4))

