function u = nonlinheat(dx,aFLAG);
% NONLINHEAT  Approximate nonlinear heat eqn 
%     u_t = e^{-5 arctan(u)} u_xx - a(x,t) u_x - 1
% using adaptive time-stepping explicit method.
% Example: >> nonlinheat(0.05,0);  % no advective term (a(x,t) = 0)
%   >> nonlinheat(0.05,1); % a(x,t) = 8 t cos(pi x)/(1+8t)

% ELB 3/28/07
tf = 0.5; 
J = floor(1/dx);  dx = 1/J;  x = 0:dx:1;  
Nmax = 500000; % stops if too many steps
t(1) = 0;  % t(n) is time at start of step n

u0 = x.*(1-x); u = u0; % length(u) = J+1
n = 0;
while ((n < Nmax) && (t(n+1) < tf))
    n = n+1;
    b = exp(-5 * atan(u(2:end-1)));
    if aFLAG == 1    
      a = 8 * t(n) * cos(pi*x(2:end-1)) / (1 + 8 * t(n));
      dt = dx/(max(2*b/dx + abs(a)));  % use (2.149)
    else
      dt = (dx^2)/(2*max(b));  % use (2.171)
    end
    t(n+1) = t(n) + min(dt, tf-t(n)); % don't go past final time
    dt = t(n+1) - t(n);   mu = dt/(dx^2);
    v = u(2:end-1);
    u(2:end-1) = v + mu * b .* (u(1:end-2) - 2*v + u(3:end)) - dt;
    if aFLAG == 1    
       for j=2:J
          if a(j-1) >= 0,  du = u(j) - u(j-1);  % upwind difference u
          else,  du = u(j+1) - u(j);  end
          u(j) = u(j) - dt * a(j-1) * du / dx;
       end
    end
end
disp(['number of steps = ' num2str(n) '; final time = ' num2str(t(n+1))])

figure(1), plot(x,u,'--',x,u0,'-'), xlabel x, ylabel u

figure(2), plot(t(1:end-1),diff(t),'.')
xlabel t, ylabel('time step \Delta t')