% type (or copy) the following lines of matlab
clear
% You might want to increase T
T=0.01;
% You might want to increase J
J=32;
h=1/J;
% looking at mu in the range (0,1) does anything interesting happen?
% Hint what is the critical mu bigger than which the behavious is strange
mu=1.0;
k=mu*h.*h;
N=T/k;
% The N above may be a non-integer, so below we round it and redefine
% the temporal discretization parameter.
N=round(N);
k=T/N;
x=[0:h:1]';
t=[0:k:T]';
U=zeros([J+1,N+1]);
U(:,1)=x.*(1-x);
for n=1:N
U([2:J],n+1)=mu.*U([1:J-1],n)+(1-2.*mu).*U([2:J],n)+mu.*U([3:J+1],n);
end
surf(x,t,U')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Forwards heat equation
%
% Used in 4H Numerical Analysis
%
% J.F.Blowey
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
