% 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 demo1

a  = 1;        % First t-value.
b  = 2.5;      % Final t-value.
y0 = 0.25;     % Initial data.
N  = 10;       % Number of discretization points.
h  = (b-a)/N;  % Step size.
f  = inline('t+y','t','y'); % Define the function "f" in the ODE.

% Draw the gradient field.
nqt  = 20; % Number of points in t (for plotting only).
nqy  = 20; % Number of points in y (for plotting only).
ymin = -1; % Lower limit of y-axis (for plotting only).
ymax = 7;  % Upper limit of y-axis (for plotting only).
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.
% We use the build-int Matlab function "ode45".
% The tolerance is set to 10^(-10).
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 "ad hoc" 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