function [x,y] = pwLadapt(f,L,fL,R,fR,tol)
% Produces column n-vectors x and y so that plot(x,y) is a good
% rendition of the function f across the interval [L,R].
% fL is the value of f at L and fR is the value of f at R.
% tol > 0 determines how close the piecewise linear interpolant
% is to f. The x and y vectors may be as long as 2(R-L)/tol so
% tol should not be too small.
if (R-L) <= tol
% The linear interpolant is acceptable...
x = [L;R]; y = [fL;fR];
else
m = (L+R)/2; fm = f(m);
d_min = Pt2LineSeg(m,fm,L,fL,R,fR);
if d_min <= tol
% The linear interpolant is acceptable...
x = [L;m;R]; y = [fL;fm;fR];
else
% Produce left and right piecewise linear interpolants...
[xL,yL] = pwLadapt(f,L,fL,m,fm,tol);
[xR,yR] = pwLadapt(f,m,fm,R,fR,tol);
% and synthesize...
x = [ xL;xR(2:length(xR))];
y = [ yL;yR(2:length(yR))];
end
end