clear all  % unassign all symbols
close all  % close all graphics windows
a = -1;
b =  1;
T = 1;
n = 20;
m = 50;
dx = (b - a) / n;
dt = T / m;
r = dt/dx^2
f = @(x) 1 - x.^2;
alpha = @(t) 0;
beta = @(t) 0;
x = a:dx:b;
t = 0:dt:T;
u = zeros(m+1, n+1);
d = ones(n+1, 1);     %   diagonal will be (1-2*r)*d
M = spdiags([ r*d  (1-2*r)*d r*d], -1:1, n+1, n+1);
M(1,1) = 0;
M(1,2) = 0;
M(end,end) = 0;
M(end,end-1) = 0;
u(1,:) = f(x);
q = zeros(n+1, 1);
for i = 1:m
    q(1) = alpha(t(i+1));
    q(end) = beta(t(i+1));
    tmp = M*u(i,:)' + q;
    u(i+1,:) = tmp';
end

surf(x,t,u);
colormap jet;
axis square;
%shading interp;
%colorbar;



