function [y,yhat] = stp(x, M)  

% [0,1] periodic step function with jumps at x = 0, 0.5.
% Returns N-vector of gridpoint values y at gridpoints x
% uniformly spaced over [0,1) and vector of Fourier coeffs yhat(M)
% ordered like the entries of M.

  y = sign(x - 0.5);   % y = -1 for 0<x<0.5, 1 for 0.5<x<1, 0 at x=0.5
  y(1) = 0;            % Set y = 0 at x=0 for symmetry with 0.5
  yhat = -2i*mod(M,2)./(pi*M);