function w = nonlinBVP(N,M);
% NONLINBVP: w=nonlinBVP(N,M)
%    Solves
%      u'' + .5 u^2 = -10*sin(pi x), u(0)=u(1)=0
%    by finite differences with Newton-Raphson iteration.
% Try "max(nonlinBVP(200,M))" for M=0,1,2,4,6 to see convergence.
% (ELB 4/26/02)

h = 1/N;  x=0+h:h:1-h; % note we want x_1=h and x_(N-1)=1-h
L=-2*eye(N-1) + diag(ones(N-2,1),1) + diag(ones(N-2,1),-1);
L=(1/h^2)*L;
B=-10; Bsinpix=B*sin(pi*x); 
w=(-1/(pi*pi))*Bsinpix; dw=zeros(size(w'));

plot([0 x 1], [0 w 0],'m:'); hold on;
title(['NONLINBVP: "First guess" dotted, soln solid. N = ' num2str(N)...
      ' intervals and ' num2str(M) ' Newton iters.']);
for l=1:M
   F=(1/h^2)*([w(2:N-1) 0] - 2 * w + [0 w(1:N-2)]) + .5*w.^2 - Bsinpix;
   DF=L+diag(w);
   dw=-DF\ F'; % Newton-Raphson step.  Note backslash.
   w=w+dw';
end
plot([0 x 1], [0 w 0],'b'); hold off;