function yp = sinsum(M)
%SINSUM  Write 1 as sum of sin series
% Try "sinsum(40);"
% (Ed Bueler; 1/18/02)
%
% Uses M terms in Fourier sum.

% create x coordinates: 501 points
xp=0:.002:1;

% sum Fourier series
yp=zeros(size(xp));
for j=M:-1:1  % count down for less rounding error
   JPI=(2*j-1)*pi;
   yp=yp+ (4/JPI)*sin(JPI*xp);
end;

% display
figure(1);
plot(xp,yp);
title('Function f(x)=1, as a sum of a sin series.');