function pts = lipgraphmesh(a,b,deg,g)

% by Federico Piazzon and Marco Vianello, University of Padova 
% December 3, 2013

% computes a sub-optimal polynomial mesh of cardinality O(deg^3) and 
% constant C close to 2 on a Lipschitz graph domain, a<=x<=b, 0<=y<=g(x), 
% with g(x) positive and Lipschitz continuous

% a call has the form pts = lipgraphmesh(a,b,deg,@f) where f.m 
% is a user defined Matlab function  

% warning: the number of points can be huge (and thus the algorithm slow), 
% for example of the size of 176*(b-a)*(L/min(g))*deg^3 
% where L is the Lipschitz constant 


% computes the minimum of g 
[mpt,m]=fminbnd(@(x) g(x),a,b);

% estimates the Lipschitz constant of g
n=100; go=1; 
h=(b-a)/(n-1);
x=linspace(a,b,n);
Lold=max(abs(g(x(2:end))-g(x(1:end-1)))/h);
while go==1
n=2*n;
h=(b-a)/(n-1);
u=linspace(a,b,n);
Lnew=max(abs(g(u(2:end))-g(u(1:end-1)))/h);
go=abs(Lnew/Lold-1)>0.001;
Lold=Lnew;
end;
L=Lnew;

% computes a sub-optimal polynomial mesh with constant 
% C = C(t) = sqrt(2)/(1-4*t*exp(2*t)), t < 0.175 
% and cardinality O(deg^3/t)
% the following value of t guarantees that C is close to 2 
t=0.0644;
ptsx=[ ];ptsy=[ ];
% computes w such that 4/w is an estimate of the Markov constant 
% of the graph domain 
w=min(m/max(2*L,sqrt(1+L^2)),(b-a)/2);
% computes the x-grid
cardx=ceil((b-a)*4/w*deg^2/(sqrt(2)*t))+1;
xx=linspace(a,b,cardx);
% computes the polynomial mesh 
ind=(0:2*deg-1)';
chebpts=cos((pi*ind+pi/2)/(2*deg));
for j=1:length(xx)-1
[z(j),M(j)]=fminbnd(@(x) -g(x),xx(j),xx(j+1));
M(j)=-M(j);
ptsy=[ptsy;M(j)/2*(1+chebpts)];
ptsx=[ptsx;repmat(z(j),2*deg,1)];
end;

pts=[ptsx ptsy];



