% August 28, 2014
% P.G. Martinsson
%
% This is a demo for APPM2360.
% It illustrates a few different methods for numerically solving the ODE
%
%    dy/dt = -t/y
%     y(a) = y0
%
% over the interval I = [a,b].

function demo2

% Set parameters:
a  = 0;
b  = 0.98;
y0 = 1;
N  = 10;
h  = (b-a)/N;
f  = inline('-t./y','t','y');

% Draw the gradient field.
nqt  = 20;
nqy  = 20;
ymin = -1.1;
ymax = 1.1;
figure(1)
[tt,yy] = meshgrid(linspace(a-0.5,b+0.5,nqt),linspace(ymin,ymax,nqy));
hold off
quiver(tt,yy,ones(nqt,nqy),f(tt,yy),'g')
title('The differential equation dy/dt = -t/y','FontSize',16)
hold on
keyboard

% Construct an approximate solution using a GOOD ode solver.
opts = odeset('AbsTol',1e-10);
[tt,yy] = ode45(f,[a,b],y0,opts);
plot(tt,yy,'k')
axis([a-0.5,b+0.5,ymin,ymax])
y_at_b = yy(end);
legend('','Reference solution')
keyboard

% Construct an approximation using forwards Euler.
tt = linspace(a,b,N+1);
yy = zeros(size(tt));
yy(1) = y0;
for n = 1:N
  yy(n+1) = yy(n) + h*f(tt(n),yy(n));
end
plot(tt,yy,'r-x')
legend('','Reference solution','Forwards Euler')
keyboard
fprintf('Error in forwards Euler             = %15.10f\n',abs(yy(end) - y_at_b))

% Construct an approximation using the 2nd order forwards scheme.
tt = linspace(a,b,N+1);
yy = zeros(size(tt));
yy(1) = y0;
yy(2) = y0 + h*f(tt(1),yy(1));
for n = 2:N
  yy(n+1) = yy(n-1) + 2*h*f(tt(n),yy(n));
end
plot(tt,yy,'b-x')
fprintf('Error in ad hoc second order scheme = %15.10f\n',abs(yy(end) - y_at_b))

% Construct an approximation using the 2nd order Runge-Kutta.
tt = linspace(a,b,N+1);
yy = zeros(size(tt));
yy(1) = y0;
for n = 1:N
  k1 = f(tt(n),yy(n));
  k2 = f(tt(n) + 0.5*h, yy(n) + 0.5*h*k1);
  yy(n+1) = yy(n) + h*k2;
end
plot(tt,yy,'m-x')
fprintf('Error in second order Runge-Kutta   = %15.10f\n',abs(yy(end) - y_at_b))

% Construct an approximation using the 4th order Runge-Kutta.
tt = linspace(a,b,N+1);
yy = zeros(size(tt));
yy(1) = y0;
for n = 1:N
  k1 = f(tt(n),yy(n));
  k2 = f(tt(n) + 0.5*h, yy(n) + 0.5*h*k1);
  k3 = f(tt(n) + 0.5*h, yy(n) + 0.5*h*k2);
  k4 = f(tt(n) + h, yy(n) + h*k3);
  yy(n+1) = yy(n) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end
plot(tt,yy,'c-x')
fprintf('Error in fourth order Runge-Kutta   = %15.10f\n',abs(yy(end) - y_at_b))
legend('','Reference solution','Forwards Euler','2nd order ad hoc','2nd order Runge-Kutta','4th order Runge-Kutta')
hold off

return