function U = thomas(J,N,tf)
% THOMAS Solves  u_t = (1+2x) u_xx  and is about as brief as possible.
dx = 1/J;  mu = (tf/N) / (dx*dx)
x = 0:dx:1;  U = x .* (1-x);
for n=1:N
   U(2:J) = U(2:J) + mu * (1+2*x(2:J)) .* (U(3:J+1) - 2 * U(2:J) + U(1:J-1));
end
plot(x,x.*(1-x), x,U,'--');