My algorithm for the heat equation is unstable

I have implemented the 2D heat equation with what I thought was the Crank-Nicolson algorithm in the following way:

function[newDomain] = heatCN(domain, dt, Nx, Ny, hx, hy, spanX, spanY)
aux = domain;
newDomain = domain;
for i=2:Nx-1
for j=2:Ny-1
derX1 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
derY1 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
aux(i,j) = derX1 + derY1 + domain(i,j);
end;
end;

for i=2:Nx-1
for j=2:Ny-1
derX2 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
derY2 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
derXT = dt * (-2.0 * aux(i,j) + aux(i+1, j) + aux(i-1,j)) / hx^2;
derYT = dt * (-2.0 * aux(i,j) + aux(i, j+1) + aux(i, j-1)) / hy^2;
newDomain(i,j) = domain(i,j) + 0.5 * (derX2 + derY2 + derXT + derYT);
end;
end;
end;

From my knowledge, the CN method should be unconditionally stable, meaning that for any (decent) time step it should give me the right output. But the behavior I see with this program is unstable.

• My advice would be to try solving this in 1D first. You should also state your boundary conditions and initial conditions. If you have initial conditions with very sharp features then this can causes oscillations which may cause divergence. Finally, you are more likely to get a good response when you post your method mathematically rather than just showing code, for example a similar 'debugging' question where the method was emphasised and not the code, scicomp.stackexchange.com/questions/5355/… – boyfarrell Dec 3 '13 at 0:36
• Thank you, @boyfarrell. The method I am trying to implement is just the standard 2D-diffusion Crank-Nicolson (en.wikipedia.org/wiki/…) I am on the other hand not using a matrix to solve the system, but trying to do this on the grid, thus saving space. I am though not sure if the idea of calculating the next time step in aux and then using that to calculate the new answer is the right one.. – Daniel Dec 3 '13 at 8:41

• @Daniel No local process can perform the solve, since the operation (the matrix inverse) is dense. You can apply this dense operation fast (in $O(n)$ time) using a hierarchical method like multigrid. If you are programming in C or C++, you may want to check out a library like PETSc, which provides many efficient solvers. – Jed Brown Dec 3 '13 at 14:54