% BOB Solves the PDE   u_t = (1+2x) u_xx

J=10;
N=1000;
tf=0.01;

dx = 1/J;
dt= tf/N;

mu = dt / (dx*dx);

x = 0:dx:1;

Uold = zeros(size(x));
for j=1:J+1
    Uold(j) = x(j) * (1 - x(j));
end
plot(x,Uold); hold on

Unew = Uold;
for n=1:N
    Uold(1) = 0;  Uold(J+1) = 0;
    for j=2:J
        Unew(j) = Uold(j) + mu * (1+2*x(j)) * (Uold(j+1) - 2 * Uold(j) + Uold(j-1));
    end
    Uold = Unew;
end

plot(x,Unew,'--'); hold off