function err=boundaryCN(nu,J,tf,BCmeth);
% BOUNDARYCN Solves heat equation problem
%    u_t = u_xx, u_x(0,t)=3 cos(3) e^{-9t}, u(1,t)=0, u(x,0)=sin(3(x-1))
% by Crank-Nicolson method.  Allows choice of boundary condition
% approximations; either (2.103) or (2.114).  Computes error at last time;
% note exact solution is u(x,t)=e^{-9t} sin(3(x-1)).
% Format:
%   err=boundaryCN(nu,J,tf,BCmeth)
% where err = max(abs( (approx. at tf)-(exact at tf) ) ),
%        nu = dt/dx^2
%         J = number of subintervals for x, on [0,1]
%        tf = final time
%    BCmeth = either '2.103' or 2.114'
%
% Examples:
%   >> boundaryCN(5,100,0.1,'2.103')
%   >> boundaryCN(5,100,0.1,'2.114')
% ELB 3/21/05

dx=1/J; x=0:dx:1;  % mesh in space
dt=nu*dx^2;  M=ceil(tf/dt);  dt=tf/M;  % mesh parameters in time

% C.-N. matrix w/o b.c. implementation; A is J by J
d=repmat(1+nu,1,J-1); % diagonal elements
sp=repmat(-0.5*nu,1,J-2); % super-diagonal elements
sb=repmat(-0.5*nu,1,J-1); % super-diagonal elements
A=sparse([2:J, 2:J-1, 2:J],[2:J, 3:J, 1:J-1],[d sp sb],J,J);

% b.c. implementation;  note: "  if  BCmeth=='2.103'  " does not work!
if strcmp(BCmeth,'2.103')
    A(1,1)=-1;  A(1,2)=1;  goodflag=false;
elseif strcmp(BCmeth,'2.114')
    A(1,1)=1+nu;   A(1,2)=-nu;   goodflag=true;
else, error('boundary method must be "2.103" or "2.114"'), end
%figure(2), spy(A), figure(1)

U=uu(x,0)';   b=zeros(J,1);  % column vectors
for n=1:M
    if goodflag  % (2.114)
        b(1)=(1-nu)*U(1)+nu*U(2)-dx*nu*( g((n+1)*dt) + g(n*dt) );
    else  % (2.103)
        b(1)=dx*g((n+1)*dt);  end
    b(2:J)=0.5*nu*U(1:J-1) + (1-nu)*U(2:J) +0.5*nu*U(3:J+1);
    U(1:J)=A\b;
end
uex=uu(x,tf)';   err=max(abs(U-uex));
plot(x,U,'.',x,uex,x,uu(x,0),'--')
legend(['approx soln at t_f=' num2str(tf)],...
    ['exact soln at t_f=' num2str(tf)],'initial condition')


function y=uu(x,t);
y=exp(-9*t)*sin(3*(x-1));

function z=g(t);
z=3*cos(3)*exp(-9*t);