% This script solves the one-dimensional BVP
%    -a u'' + b u' = 0, 0 < x < 1
%   u(0) = 0, u(1) = 1
% using central differences for the diffusion term and
% central differences, forward differences and backward
% differences for the convection term
%
% By: Martijn Anthonissen, m.j.h.anthonissen@tue.nl

a = 1; b = 10;

n = 20;
h = 1/(n+1);
x = [1:n]*h;

% Central differences for both derivatives
e = ones(n,1);
Adiff = 1/h^2 * spdiags([-e 2*e -e], -1:1, n, n);
AconvCD = 1/(2*h) * spdiags([-e 0*e e], -1:1, n, n);
A = a * Adiff + b * AconvCD;
f = zeros(n, 1);
f(n) = a/h^2 - b/(2*h);

uexact = (exp(b/a * x) - 1) / (exp(b/a) - 1);
u = A\f;

plot([0 x 1], [0 uexact 1], '-')
hold on
plot(x, u, 'o')
title('Central differences')

% Forward differences
AconvFD = 1/h * spdiags([0*e -e e], -1:1, n, n);
A = a * Adiff + b * AconvFD;
f = zeros(n, 1);
f(n) = a/h^2 - b/h;

u = A\f;

figure
plot([0 x 1], [0 uexact 1], '-')
hold on
plot(x, u, 'o')
title('Forward differences')

% Backward differences
AconvBD = 1/h * spdiags([-e e 0*e], -1:1, n, n);
A = a * Adiff + b * AconvBD;
f = zeros(n, 1);
f(n) = a/h^2 + b/h;

u = A\f;

figure
plot([0 x 1], [0 uexact 1], '-')
hold on
plot(x, u, 'o')
title('Backward differences')