[Edit: Correct "x(i)=pos(1) & y(j)=pos(2)" to "x(i)==pos(1) & y(j)==pos(2)".]
[Edit: Fixed the update formula.]
Am I correct to infer from your code that there is a heat source at location pos?
Am I correct that u(i,j) is the teperature at position x(i),y(j)?
What is the temperature at the source? I do not see that in the code.
Am I correct that temperature is zero along the edges?
I do not udenrstand why you use the exponential function in your solution. I realize that the exponential is the analytic solution. IT seems you are also partly using the finite difference method for the second derivative, but I don;t understand how you are combining these approaches.
I would use a simple forward finite difference method with a small time step and see what happens. I realize this is not as elegant as Crank-Nicholson or other higher order approaches, but it is simple and should work, as long as the time step is small enough compared to the spatial grid.
My loops would be something like this:
uold=zeros(nx,ny); unew=uold;
for n=1:nt
for i=2:nx-1
for j=2:ny-1
if x(i)==pos(1) & y(j)==pos(2)
unew=Tsource;
else
unew(i,j)=uold(i,j)+(uold(i-1,j)+uold(i+1,j)+uold(i,j-1)+uold(i,j+1)-4*uold(i,j))*k*dt/h^2;
end
end
end
uold=unew;
%plot unew, and wait a bit before continuing
end
The loop above does not update the edge elements, so they will stay at u=0 (their initial state) throughout.
Good luck.
