% script: exer016.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 = 1;
fun = 'fcn016';
numSteps = 30*[1,2,4,8];
hSteps = (tLast - t0)./numSteps;
lenSteps = length(hSteps);
lerrors = zeros(lenSteps,3);
% store stepsizes and errors in this matrix
errors(:,1) = hSteps;
% main  -- we use the infinity norm for error
for jj = 1:lenSteps
  tInterval = [t0:hSteps(jj):tLast];
  ySoln = soln016(tInterval);
  yApprox = Euler(fun,t0,y0,hSteps(jj),numSteps(jj));
  errors(jj,2) = norm(yApprox - ySoln,inf);
  yApprox = ERK4(fun,t0,y0,hSteps(jj),numSteps(jj));
  errors(jj,3) = norm(yApprox - ySoln,inf);
end
disp('Table with columns [time stepsize, Euler error, ERK4 error]:');
format long
disp(errors)
disp('Table with columns of rates of decrease of error [stepsize, Euler, ERK4]');
disp(errors(1:lenSteps-1, 1:3)./errors(2:lenSteps, 1:3))
format

