function [t,err]=nonlinverif(dx,tf)
% NONLINVERIF  Function to approximate solution of nonlinear heat
% equation
%   u_t = e^{-30 arctan(u)} u_xx + f(x,t),  u(0,t)=u(1,t)=0
% using adaptive time-stepping explicit method.  Compares to exact soln
%   u(x,t) = e^{-t/3} sin(pi x).
% (In the case f(x,t) is computed from exact solution and PDE; i.e.
% it is manufactured.)  Plots and displays error at tf in figure(1).
% Example:
%   >> [t,err]=nonlinverif(0.02,5.0);
%   >> figure, plot(t(1:end-1),t(2:end)-t(1:end-1),'.')
%   >> xlabel t, ylabel('time step \Delta t')
%   >> figure, plot(t,err), xlabel t, ylabel('error')
% ELB 3/30/05

Nmax=500000; % stops if too many steps
t(1)=0; % t(n) = time at end of (n-1)st step
err(1)=0; % err(n) = error at end of (n-1)st step
x=0.0:dx:1.0;
[ignor,u]=fu(x,0);  % get initial condition from exact
for n=1:Nmax
    a=exp(-30*atan(u(2:end-1)));
    dt=(dx^2)/(2*max(a));  % from (2.171)
    if t(n)+dt > tf  % if step unnecessarily long at end
        dt=tf-t(n); end
    nu=dt/(dx^2);    t(n+1)=t(n)+dt;    v=u(2:end-1);
    [f,uex]=fu(x,t(n));  % compute exact soln and source term
    err(n+1)=max(abs(u-uex)); % compare to exact
    u(2:end-1)=v+nu*a.*(u(1:end-2)-2*v+u(3:end)) + dt*f(2:end-1);
    if t(n)+dt >= tf, break, end  % if done
end
if n>=Nmax, warning('max number of iterations reached'), end
disp(['total number of steps: ' num2str(n+1)])
% plot exact and approximate at final time
[ignor,uex]=fu(x,tf);
figure(1), plot(x,u,'.',x,uex), xlabel x, ylabel u
title(['time = ' num2str(tf) '; err = ' num2str(max(abs(u-uex)))])
legend('approximate','exact')
        
% compute source term and exact soln together; manufactured soln
function [f u]=fu(x,t);
u=exp(-t/3)*sin(pi*x);
f=( -(1/3) + pi^2 * exp(-30 *atan(u)) ).* u;