% Program 8.2 Backward difference method for heat equation
% input: space interval [xl,xr], time interval [yb,yt],
% number of space steps M, number of time steps N
% output: solution w
% Example usage: w=heatbd(0,1,0,1,0.02,0.02)
function w=heatbd(xl,xr,yb,yt,h,k)
f = @(x) 2*cosh(x); % Oppgave 1a
l = @(t) 2*exp(2*t);
r = @(t) (exp(2)+1)*exp(2*t-1);
D=2; % diffusion coefficient
M=(xr-xl)/h; % number of space steps
n=(yt-yb)/k; % number of time steps
m=M-1;
sigma=D*k/(h*h);
a=diag(1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);
a=a+diag(-sigma*ones(m-1,1),-1); % define matrix a
lside=l(yb+(0:n)*k); rside=r(yb+(0:n)*k);
w(:,1)=f(xl+(1:m)*h)'; % initial conditions
for j=1:n
w(:,j+1)=a\(w(:,j)+sigma*[lside(j);zeros(m-2,1);rside(j)]);
end
w=[lside;w;rside]; % attach boundary conds
x=(0:m+1)*h;t=(0:n)*k;
mesh(x,t,w') % 3-D plot of solution w
view(60,30);axis([xl xr yb yt -1 60])