% quarter wave sine series for solution to
%-u''= f(x)=x on [0,1], u(0)=0, du/dx(1)=0
n=input('Choose N');
clear c;
for j=1:n
    c(j)=32*(-1)^(j+1)/(((2*j-1)^4)*pi^4);
end;
xgrid=0:0.01:1;
hold off
fn=zeros(size(xgrid));
for j=1:n
    fn=fn+ c(j)*sin((2*j-1)*pi*xgrid/2);
end
plot(xgrid,-xgrid.^3/6+xgrid/2, xgrid, fn);