function x = meshgen_mmpde5(N,y,rho,xi,dt,tau)
%MESHGEN_MMPDE5 generates a new mesh of N points using the MMPDE5 algorithm.
%
% X = MESHGEN_MMPDE5(N,Y,RHO) generates a new mesh of N points using
% the MMPDE5 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_MMPDE5(N,Y,RHO,XI)
% X = MESHGEN_MMPDE5(N,Y,RHO,XI,DT)
% X = MESHGEN_MMPDE5(N,Y,RHO,XI,DT,TAU)
%
% The optional arguments are XI, DT and TAU. Their default values are
% XI(i) = (i-1)/(N-1), i = 1,...,N,
% DT = 0.1, and TAU = 1.
%
% 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);
A = sparse(N,N);

% 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
   dt = 0.1;
end
if nargin < 6
   tau = 1.0;
end

% Assemble the coefficient matrix and the right-hand side term 

dt1 = dt/tau;
A(1,1) = 1; 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));
   A(i,i-1) = -0.5*(rho(i)+rho(i-1))*dxi*dxi1*dt1/rho(i);
   A(i,i+1) = -0.5*(rho(i)+rho(i+1))*dxi*dxi2*dt1/rho(i);
   A(i,i) = 1.0 - A(i,i-1) - A(i,i+1);
   x(i) = y(i);
end
A(N,N) = 1; x(N) = y(N);

% solve the system for the new mesh

x = A\x;
x(1) = y(1); x(N) = y(N);

% end of meshgen_mmpde5
