function x = meshgen_mmpde5ex(N,y,rho,xi,dt,tau)
%MESHGEN_MMPDE5EX generates a new mesh of N points using the MMPDE5
% (Explicit) algorithm.
%
% X = MESHGEN_MMPDE5EX(N,Y,RHO) generates a new mesh of N points using
% the MMPDE4 algorithm. Y is the old mesh of N points. RHO is the adaptation
% function defined on the old mesh Y. 
%
% Usage with optional arguments:
% X = MESHGEN_MMPDE5EX(N,Y,RHO,XI)
% X = MESHGEN_MMPDE5EX(N,Y,RHO,XI,DT)
% X = MESHGEN_MMPDE5EX(N,Y,RHO,XI,DT,TAU)
%
% The optional arguments are XI and DT. Their default values are
% XI(i) = (i-1)/(N-1), i = 1,...,N and DT is calculated according to
% stability consideration.
%
% XI specifies the computational mesh of N points.
%
% DT sets the pseudo-time step size.
%
% TAU sets the value of parameter tau.
%

%
% 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.
%

% declare arrays

x = zeros(N,1);

% Initialization

if nargin < 3
   error('The first three arguments must be given.')
end
if nargin < 4
   xi = linspace(0,1,N)';
end
if nargin < 5
   temp = 1.0;
   for i=2:N-1
      dxi = 2/(xi(i+1)-xi(i-1));
      dxi1 = 1/(xi(i)-xi(i-1));
      dxi2 = 1/(xi(i+1)-xi(i));
      temp1 = 0.5*(rho(i)+rho(i-1))*dxi1;
      temp2 = 0.5*(rho(i)+rho(i+1))*dxi2;
      temp1 = dxi*(temp1+temp2)/rho(i); 
      temp1 = 1.0/temp1;
      if temp1 < temp, temp = temp1; end
   end
   dt = temp*0.8;
end
if nargin < 6
   tau = 1.0;
end

% explicit euler integration

dt1 = dt/tau;
x(1) = y(1);
for i = 2:N-1
   dxi = 2/(xi(i+1)-xi(i-1));
   dxi1 = 1/(xi(i)-xi(i-1));
   dxi2 = 1/(xi(i+1)-xi(i));
   temp1 = 0.5*(rho(i)+rho(i-1))*(y(i)-y(i-1))*dxi1;
   temp2 = 0.5*(rho(i)+rho(i+1))*(y(i+1)-y(i))*dxi2;
   x(i) = y(i)+dt1*dxi*(temp2-temp1)/rho(i);
end
x(N) = y(N);

% end of meshgen_mmpde5ex
