function [U,err,time]=molmanuf(J,tf,method)
% MOLMANUF  Function which demonstrates method of lines on a nonlinear
% heat equation with a manufactured solution:
%     u_t = 0.1 u_xx - 2 u^2 + f(x,t),   (x,t) in (0,1) x (0,tf),
% where
%     f(x,t) = 0.1 (9 pi^2/4 - 1) e^{-0.1 t} cos(3 pi x/2)
%              + 2 e^{-0.2 t} cos^2(3 pi x/2),
% with boundary conditions  u_x(0,t)=0,  u(1,t)=0   and initial condition
%     u(x,0) = cos(3 pi x/2).
% Uses various built-in ODE solvers to solve semidiscretized problem. 
% Note exact solution is    u(x,t) = e^{-0.1 t} cos(3 pi x/2).
% Form:
%    [U,err,time]=molmanuf(J,tn,method)
% where:
%        U = approximate solution at times in tn
%            ( row k of U is length J+1; U(k,:) is approx soln at tn(k) )
%      err = maximum error compared to exact solution at t=tf
%     time = time in seconds for ode solve
%        J = number of subintervals
%       tf = final time
%  [method = 'ode23s' (default) or 'ode45' or 'ode23sSMALLTOL'; OPTIONAL;
%            Note: default RelTol=1e-3, AbsTol=1e-6  for both ode23s and
%            ode45.  For 'ode23sSMALLTOL', RelTol=1e-6, AbsTol=1e-9
% Example 1; shows error optimal given spatial discretization
%    >> J=[10 20 30 40 60 100 140 200 250 300 400 500];  tf=0.1;
%    >> for k=1:12, [U,err(k),tm(k)]=molmanuf(J(k),tf); tm(k), end
%    >> loglog(1./J,err,'o'), xlabel('\Delta x'), ylabel('error')
%    >> p=polyfit(log(1./J),log(err),1)
% Example 2; nonstiff solver takes longer w. no increase in accuracy:
%    >> [U,er,tm]=molmanuf(700,0.1,'ode23s');  er, tm  % a few seconds
%    >> [U,er,tm]=molmanuf(700,0.1,'ode45');  er, tm  % a few minutes
% Example 3; takes several minutes, but one can see effect of tolerances:
%    >> [U,er,tm]=molmanuf(2000,0.1,'ode23s');  er, tm
%    >> [U,er,tm]=molmanuf(2000,0.1,'ode23sSMALLTOL');  er, tm
% ELB 4/11/05
% See also: MOLEXAMPLE

dx=1/J;   x=0:dx:1-dx;   U0=cos(3*pi*x/2);
if nargin==2
    tic, [t,U]=ode23s(@molderivs,[0.0 tf],U0); time=toc;
elseif nargin==3
    if strcmp(method,'ode45')
        tic, [t,U]=ode45(@molderivs,[0.0 tf],U0); time=toc;
    elseif strcmp(method,'ode23s')
        tic, [t,U]=ode23s(@molderivs,[0.0 tf],U0); time=toc;
    elseif strcmp(method,'ode23sSMALLTOL')
        options=odeset('RelTol',1e-6,'AbsTol',1e-9);
        tic, [t,U]=ode23s(@molderivs,[0.0 tf],U0,options); time=toc;
    else, error('method must be "ode45", "ode23s", "ode23sSMALLTOL"'), end
else, error('requires 2 or 3 input arguments'), end
err=max(abs( U(end,:) - exp(-0.1*tf)*U0 ));

function dUdt=molderivs(t,U);
% U,dUdt are a column vectors;  t not used
J=length(U);  dx=1/J;  x=0:dx:1-dx;
nu=0.1/(dx*dx);  dUdt=zeros(J,1);
u=exp(-0.1*t)*cos(3*pi*x/2);  % exact solution; only used to build f
f=0.1*(9*pi^2/4 - 1)*u + 2*u.^2;
dUdt(1)=-2*nu*U(1)+2*nu*U(2)-2*U(1)^2 + f(1);
for k=2:J-1
    dUdt(k)=nu*U(k-1)-2*nu*U(k)+nu*U(k+1)-2*U(k)^2 + f(k);   end
dUdt(J)=nu*U(J-1)-2*nu*U(J)-2*U(J)^2 + f(J);