  function [x,y] = pwLStatic(fname,M2,alpha,beta,delta)
%
% Pre:
%  fname        string that names an available function f(x).
%               Assume that f can take vector arguments.
%  M2           an upper bound for|f"(x)| on [alpha,beta].
%  alpha,beta   scalars with alpha<beta
%  delta        positive scalar
%
% Post:
%  x,y          column n-vectors with the property that y(i) = f(x(i)),
%               i=1:n. The piecewise linear interpolant L(x) of this
%               data satisfies |L(x) - f(x)| <= delta, where
%               x = linspace(alpha,beta,n).
%
  n = max(2,ceil(1+(beta-alpha)*sqrt(M2/(8*delta))));
  x = linspace(alpha,beta,n)';
  y = feval(fname,x);
  
