% script: exerVarTest.m
% description: solves ODEs with various methods

 % programmer: t. shores
% last rev: 4/7/16

global gError gxCount;

% Octave and Matlab compatible file
disp('Exercise 1, NumAnal Seminar 2016:')
disp('dy/dt = -t*y*y,  0<=t<=6')
disp('y(0) = 1')
%disp('Curtiss-Hirschfeld:')
%disp('dy/dt = -50*(1 - cos(t)*y  0<=t<=6')
%disp('y(0) = 1')
% initialization code
t0 = 0;
tLast = 6; 
y0 = exerVarSoln(t0);
fun = 'exerVarRHS';
numSteps = 100*[1, 2, 4, 8];
hSteps = (tLast - t0)./numSteps;
lenSteps = length(hSteps);
errors = zeros(lenSteps, 4);
fcnCounts = zeros(lenSteps,3);
% store stepsizes and errors in this matrix
errors(:,1) = hSteps;
% main  -- we use the infinity norm for error
ySoln = exerVarSoln(tLast); % for end point error
for jj = 1:lenSteps
%  ySoln = exerVarSoln(linspace(t0, tLast, numSteps(jj) + 1));  % for 2-norm average error
  gxCount = 0;
%  yApprox = Euler(fun, t0, y0,hSteps(jj), numSteps(jj), 0); % for 2-norm average error
%  errors(jj, 2) = norm(yApprox - ySoln)*sqrt(hSteps(jj))/(tLast - t0); % ditto
  yApprox = Euler(fun, t0, y0,hSteps(jj), numSteps(jj));% for end point error
  errors(jj, 2) = norm(yApprox - ySoln, inf); % ditto
  fcnCount(jj, 1) = gxCount;
  gxCount = 0;
%  yApprox = ERK4(fun, t0, y0, hSteps(jj), numSteps(jj), 0); % for 2-norm average error
%  errors(jj, 3) = norm(yApprox - ySoln)*sqrt(hSteps(jj))/(tLast - t0); % ditto
  yApprox = ERK4(fun, t0, y0, hSteps(jj), numSteps(jj)); % for end point error
  errors(jj, 3) = norm(yApprox - ySoln, inf); % ditto
  fcnCount(jj, 2) = gxCount;
  gxCount = 0;
%  yApprox = ERK23(fun, t0, y0, hSteps(jj), numSteps(jj), 0); % for 2-norm average error
%  errors(jj, 4) = norm(yApprox - ySoln)*sqrt(hSteps(jj))/(tLast - t0); % ditto
  yApprox = ERK23(fun, t0, y0, hSteps(jj), numSteps(jj)); % for end point error
  errors(jj, 4) = norm(yApprox - ySoln,inf); % ditto
  fcnCount(jj, 3) = gxCount;
end
format short e
disp('Table with columns of time steps and errors:')
disp(' [Time stepsize,   Euler         ERK4,    ERK23]');
disp(errors)
disp(' ')
disp('Table with columns of rates of decrease of error:');
disp(' [Stepsize      Euler      ERK4       ERK23]');
disp(errors(1:lenSteps-1, 1:4)./errors(2:lenSteps, 1:4))
disp(' ')
format
disp('Table with function evaluations:');
disp('[Euler     ERK4     ERK23]');
disp(fcnCount)
format


