function sec2_5ex3(jmax)
%
% driver for example 2.5.3:
% generate an adaptive mesh for a given function.
%
% jmax:     the total number of mesh points. 
% METHOD = 1: for De Boor's algorithm,
%          2: for TPBVP algorithm,
%          3: for MMPDE5 algorithm,
%          4: for MMPDE5EX algorithm.
%          5: for MMPDE5XI algorithm.
%

%
% Copyright (C) 2010 Weizhang Huang and Robert D. Russell
% all rights reserved.
%
% This program is provided "as is", without warranty of any kind.
% Permission is hereby granted, free of charge, to use this program
% for personal, research, and education purposes. Distribution or use
% of this program for any commercial purpose is permissible
% only by direct arrangement with the copyright owner.
%

cpu0=clock;

   % Initialization

   % monitor = 0 (fixed mesh), 1, 2, 3, 4, 5
   % alpha_type = 1, 2, or 3: (global def, global def with flooring, or alpha = constant)
    
   monitor = 3;

   method = 5;
   tol = 1e-8;
   npde = 2;
   dt = 0.01;
   t = 0.0;
   
   % for adaptation function
   lsfitting='yes';
   alpha_type = 1;
   alpha = 1.0;

   x = zeros(jmax,1);
   xi = zeros(jmax,1);
   v = zeros(npde,jmax);
   
   x = linspace(0,1,jmax)';
   xi = linspace(0,1,jmax)';
   for j=1:jmax, v(:,j)=exactU(0.0,x(j)); end
   rho = mesh_density_fnct(npde,jmax,x,v,monitor,lsfitting,alpha_type,alpha);
   for j=1:0 % no smoothing
      rho1=rho;
	  rho = mesh_density_fnct_smoothing(jmax,rho1);
   end
   errH1_0 = interpolation_error(@exactU,jmax,npde,0.0,x,v,2,1,1,@exactDU);
   errL2_0 = interpolation_error(@exactU,jmax,npde,0.0,x,v,2,1,0,[]);

   % Main loop for convergent adaptive mesh

   for k=1:10000
      x1 = x;
      switch method
         case 1
             x = meshgen_deboor(jmax, x1, rho);
         case 2
             x = meshgen_tpbvp(jmax, x1, rho);
         case 3
             x = meshgen_mmpde5(jmax, x1, rho, xi, dt);
         case 4
             x = meshgen_mmpde5ex(jmax, x1, rho, xi, dt);
         case 5
             x = meshgen_mmpde5xi(jmax, x1, rho, xi, dt);
         otherwise
             x = meshgen_deboor(jmax, x1, rho);
      end
      for j=1:jmax, v(:,j)=exactU(0.0,x(j)); end
      rho = mesh_density_fnct(npde,jmax,x,v,monitor,lsfitting,alpha_type,alpha);
      for j=1:0 % no smoothing
         rho1=rho;
	     rho = mesh_density_fnct_smoothing(jmax,rho1);
      end

      % convergence test and info ouput

      qual = qmesh(jmax, x, rho);
      diff = norm(x1-x,inf);
      fprintf('%4d %7.3f %12.5e\n', k, qual, diff);
      Nn(k) = k;
      Qual(k) = qual;
      Diff(k) = diff;

      if diff < tol, break; end
   end
   fprintf('\n cpu time used = %e \n', etime(clock,cpu0));

   errH1 = interpolation_error(@exactU,jmax,npde,0.0,x,v,2,1,1,@exactDU);
   errL2 = interpolation_error(@exactU,jmax,npde,0.0,x,v,2,1,0,[]);

   dx=1.0;
   for j=2:jmax
      temp=x(j)-x(j-1);
	  if temp<dx, dx=temp; end
   end

   fprintf('\n %4d %12.5e %12.5e %12.5e\n',jmax,errL2,errH1,dx);
   
% plot adaptation function, mesh density, and convergence history

figure(1);
plot(x,v(1,:)','-o',x,v(2,:),'-x');
%plot(x,v(1,:)','-o');
title('(a)')
xlabel('x'); ylabel('u');
print -deps a.eps;

figure(2);
clear x1;
for i=2:jmax
   x1(i-1) = 1/(x(i)-x(i-1));
   bb(i-1) = 0.5*(x(i)+x(i-1));
end
plot(bb,x1,'-');
title('(b)')
xlabel('x'); ylabel('Mesh Density');
print -deps b.eps;

figure(3);
semilogy(Nn, Diff, '-+');
title('(c)')
xlabel('Iteration Number (k)'); ylabel('Diff Between Consecutive Meshes');
print -deps c.eps;

figure(4);
plot(Nn, Qual,'-+');
title('(d)')
xlabel('Iteration Number (k)'); ylabel('Adpt Quality Measure');
print -deps d.eps;

% -----------------------------------------------------

function u=exactU(t,x)
   u=zeros(2,1);
   u(1) = tanh(100*(x-0.3));
   u(2) = -tanh(150*(x-0.6));
      
function du=exactDU(t,x)
   du=zeros(2,1);
   du(1) = 100*(1-tanh(100*(x-0.3))*tanh(100*(x-0.3)));
   du(2) = -150*(1-tanh(150*(x-0.6))*tanh(150*(x-0.6)));
