Matlab solution for non-homogenous heat equation using finite differences

Question

Given the following PDE (non-homogenous heat equation):

$$ u_t(x,t) = c^2u_{xx}(x,t) + f(x,t) $$ $$ u(0,t) = u(l,t) = 0 $$ $$ u(x,0) = g(x) $$ $$ 0 < x < l ; t > 0 ; c > 0 $$

I wrote the following code in Matlab, to solve the problem using finite differences:

syms xj tk

% Manually define this values
c        = 9;
f(xj,tk) = xj;
g(xj)    = 0*xj;
l        = 1;
Tmax     = 0.1;

% Grid definition
Nx     = 50;          
Nt     = 50;
hx     = 1/Nx;
ht     = 1/Nt;
x      = 0:hx:l;
t      = 0:ht:Tmax;
lambda = c^2*ht/hx^2;

% Our target
u = zeros(Nx+1,Nt+1);

% Initial values
for j=1:Nx,
    u(j,1) = g(x(j));     % u(x,0) = g(x)
end

for k=1:Nx,
    u(1,k+1) = 0;         % border condition u(0,t) = 0
    for j=2:Nt,
        u(j,k+1) = u(j,k) + lambda*(u(j+1,k)-2*u(j,k)+u(j-1,k)) + ht*f(j,k); % the formula here is ok
    end
    u(Nt,k+1) = 0;        % border condition u(l,t) = 0
end

contour3(u)

For some reason that I cant't figure out, data is only appearing in the last columns and in a very strange way.

I'm guessing the implementation of the BC's are doing something nasty. But I don't see it.

Is there something that I'm missing?

Thanks in advance!

Answer 1

I doubt the last for loop should be running from $k = 1$ to $k = Nt$, and the last line in the loop should be $u(Nx, k+1) = 0$.