clear all	% unassign all symbols
close all	% close all figures
f = @(x) 1 - x.^2;
alpha = @(t) 0;
beta = @(t) sin(t);
a = -1;
b = 1;
T = pi;
n = 20;
m = 60;
dx = (b-a)/n;
dt = T/m;
r = dt/dx^2
u = zeros(m+1, n+1);
x = a:dx:b;
t = 0:dt:T;
e = ones(n+1, 1);
M = spdiags([-r*e (1+2*r)*e -r*e], -1:1, n+1, n+1);
M(1,1) = 1;
M(1,2) = 0;
M(end,end) = 1;
M(end,end-1) = 0;
u(1,:) = f(x);
for i = 1:m
    rhs = [ alpha(t(i+1)) u(i,2:end-1) beta(t(i+1)) ]';
    sol = M \ rhs;
    u(i+1,:) = sol';
end

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