function sec2_2ex1ex2(N,method)
%
% function sec2_2ex1ex2(N,method)
% It generates a convergent adaptive mesh for a given adaptation function.
% 
% choice for parameters:
% N: use a larger value if the code does not converge.
% method = 1 (for ex1) and 2 (for ex2).
%
% N:     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.
%

% Initialization

tol = 1e-8;
dt = 1.0; % dt is not used in this example.

switch method
   case 1 
      fprintf('\n\n *** De Boor''s algorithm is used *** \n\n');
   case 2 
      fprintf('\n\n *** TPBVP algorithm is used *** \n\n');
   case 3 
      fprintf('\n\n *** MMPDE5 algorithm is used *** \n\n');
   case 4 
      fprintf('\n\n *** MMPDE5EX algorithm is used *** \n\n');
   case 5 
      fprintf('\n\n *** MMPDE5XI algorithm is used *** \n\n');
   otherwise
      fprintf('\n\n *** De Boor''s algorithm is used *** \n\n');
end

x = linspace(0,1,N)';
xi = linspace(0,1,N)';
for i=1:N, rho(i) = funct_rho(x(i)); end

% Main loop

for k=1:100000
   x1 = x;
      switch method
      case 1
          x = meshgen_deboor(N, x1, rho);
      case 2
          x = meshgen_tpbvp(N, x1, rho);
      case 3
          x = meshgen_mmpde5(N, x1, rho, xi, dt);
      case 4
          x = meshgen_mmpde5ex(N, x1, rho, xi, dt);
      case 5
          x = meshgen_mmpde5xi(N, x1, rho, xi, dt);
      otherwise
          x = meshgen_deboor(N, x1, rho);
   end

   for i=1:N, rho(i) = funct_rho(x(i)); end
   % convergence test and info ouput
   qual = qmesh(N, 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
x(1) = 0;
x(N) = 1.0;

% plot adaptation function, mesh density, and convergence history

figure(1);
plot(x,rho,'-o');
title('(a)')
xlabel('x'); ylabel('Adaptation Function');
print -deps a.eps;

figure(2);
clear x1;
for i=2:N
   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;

% compute the exact equidistributing mesh and the error
% in the computed mesh
 
xi = linspace(0,1,N)';
x1 = fsolve(@(x) exact_inverse_map(x, xi), x);
err = norm(x1-x,inf);
err_cell_size = 0.0;
for i=2:(N-1)
   temp = abs(x1(i)-x(i))/(x1(i+1)-x1(i));
   if temp>err_cell_size, err_cell_size = temp; end
   temp = abs(x1(i)-x(i))/(x1(i)-x1(i-1));
   if temp>err_cell_size, err_cell_size = temp; end
end
fprintf('N = %4d err = %e err_cell_size = %e\n', N, err, err_cell_size);

% end of main 

%
% local functions defined below
%

function rho = funct_rho(x)
% this defines the adaptation function rho.
%
rho = 1+30*(1-tanh(30*(x-0.5))^2)+10*(1-tanh(10*(x-0.75))^2) ...
       +20*(1-tanh(20*(x-0.25))^2);
%
function f = exact_inverse_map(x, xi)
% this defines the exact inverse map.
sigma = 1+tanh(20*(1-0.25))+tanh(30*(1-0.5))+tanh(10*(1-0.75)) ...
        -tanh(20*(0-0.25))-tanh(30*(0-0.5))-tanh(10*(0-0.75));
f = x+tanh(20*(x-0.25))+tanh(30*(x-0.5))+tanh(10*(x-0.75)) ...
        -tanh(20*(0-0.25))-tanh(30*(0-0.5))-tanh(10*(0-0.75))-xi*sigma;
