function x = meshgen_tpbvp(N,y,rho,xi)
%MESHGEN_TPBVP generates a new mesh of N points using the TPBVP algorithm.
%
% X = MESHGEN_TPBVP(N,Y,RHO) generates a new mesh of N points using
% the TPBVP algorithm. Y is the old mesh of N points. RHO is the adaptation
% function defined on the old mesh Y. 
%
% Usage with optional argument:
% X = MESHGEN_TPBVP(N,Y,RHO,XI)
%
% The optional argument is XI. Its default value is XI(i) = (i-1)/(N-1),
% i = 1,...,N. It specifies the computational mesh of N points.
%

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

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

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;
   A(i,i+1) = 0.5*(rho(i)+rho(i+1))*dxi*dxi2;
   A(i,i) = - A(i,i-1) - A(i,i+1);
   x(i) = 0;
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_tpbvp
