function [y,kk]=p_const(x,bpts,vals)
%Computes a piecewise constant function with breakpoints BPTS 
%and corresponding values in VALS
if nargin==1,
   bpts=[min(x),(min(x)+max(x))/2,max(x)];
   vals=[-1,0,1];
end
n=length(bpts);
y=zeros(size(x));
kk=[];
for ii=1:(n-1),
   I=find(bpts(ii)<=x & x<bpts(ii+1));
   y(I)=polyval(vals(ii),x(I));
   kk(ii)=I(1);
end
J=find(bpts(n)<=x);
y(J)=polyval(vals(n),x(J));
kk(n)=J(1);
if nargout<1,
   y(kk)=NaN;
   plot(x,y,bpts,vals,'rd'), axis([min(x),max(x), 1.2*min(vals),...
         1.2*max(vals)])
end