Crank-Nicholson solution of 1D heat equation

• MATLAB
Homework Helper
I wish to numerically compute solutions of the 1D heat equation using the Crank-Nicholson scheme:

The equation is:
$$\partial_{t}u=\partial^{2}_{x}u$$

I use the discretisation:
$$u_{i+1,j}-u_{i,j}=s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})+s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})$$

Where $s=\delta t/2\delta x$. This can be rearranged into the following:
$$su_{i+1,j+1}-(2s+1)u_{i+1,j}+su_{i+1,j-1}=-su_{i,j+1}+(2s-1)u_{i,j}-su_{i,j-1}$$

This is a matrix equation which is tridiagonal and can be solved very easily using matlab's many inbuilt functions.

So suppose I wish to use the following boundary conditions $\partial_{x}u(t,0)=\partial_{x}u(t,L)=0$ how would I go about doing it? I've heard about people going from an $n\times n$ matrix to an $(n-2)\times (n-2)$, solving that and just adding in the boundary conditions after that? I tried that method but my solution blew up annoyingly.

Any suggestions on what I am doing wrong? My code is here:

Code:
%This program is meant to test out the Crank-Nicolson scheme using a simple
%nonlinear diffusion scheme.
n=2000;m=30;
t=linspace(0,20,n);% set time and distance scales
x=linspace(0,1,m);
dx=x(2)-x(1);dt=t(2)-t(1); %Increments
s=dt/(2*dx^2);%Useful for the solution.
u=zeros(n,m); %set up solution matrix
p=s*ones(1,m-1); q=-(1+2*s)*ones(1,m);
Z=diag(p,1)+diag(p,-1)+diag(q);
A=Z(2:m-1,2:m-1);
%Add in intial condition:
u(1,:)=exp(-5*(x-0.5).^2);
v=zeros(m-2,1);
%Solve the system
for i=2:n-1
%Construct the RHS for the solution
for j=2:m-1
v(j-1)=s*u(i-1,j+1)-(2*s+1)*u(i-1,j)-s*u(i-1,j-1);
end
%Solve for the new time step
w=A\v;
u(i,2:m-1)=w;
u(i,1)=u(i,2);
u(i,end)=u(i,end-1);
end

A common fall is to use first order finite difference for the derivative boundary condition. Use a 3 points formula instead. I don't know how that can be embedded in the solution for the transient problem though.

Homework Helper
It shouldn't matter with the transient problem as we're just dealing with boundary conditions. so you're suggesting I use something like:

$$\frac{\partial u}{\partial x}\approx\frac{u_{i,j+1}-u_{i,j-1}}{2\delta x}$$

That's the accuracy you're referring to right?

Homework Helper
The code is now working and I have also solved the following extension of the simple heat equation:

$$\rho c\frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\left(k(x)\frac{\partial T}{\partial x}\right)+\dot{Q}$$