% Usage: 
% y = newtfd(f,fy,fz,a,b,ya,yb,n) or [y x] = newtfd(f,fy,fz,a,b,ya,yb,n)
% 
% Finite difference method for nonlinear boundary value problems
% 
% y" = f(x,y,y')
% 
% Uses Newton's iteration for solving the nonlinear system
%
% Input:
% f - Matlab inline function f(x,y,z)
% fy - derivative of f w.r.t. y
% fz - derivative of f w.r.t. z
% a,b - interval
% ya, yb - (Dirichlet) boundary conditions
% n - number of panels
%
% Output:
% y - computed solution
% x - mesh points
%
% Examples:
% f = inline('y-y*y','x','y','z');
% fy = inline('1-2*y','x','y','z');
% fz = inline('0','x','y','z');
% [y x]=newtfd(f,fy,fz,0,1,1,4,10);
%
% f = inline('z+cos(y)','x','y','z');
% fy = inline('-sin(y)','x','y','z');
% fz = inline('1','x','y','z');
% [y x]=newtfd(f,fy,fz,0,pi,0,1,10);

function [w x] = newtfd(f,fy,fz,a,b,ya,yb,n)
h = (b - a) / n;
hh = h*h;
h2 = h / 2;
hinv2 = 0.5 / h;
% The mesh points and initial guess
x = linspace(a,b,n+1);
w = linspace(ya,yb,n+1)';
% The first and last rows of the Jacobian matrix DF(w) and the vector F(w)
J = sparse(n+1,n+1);
J(1,1) = 1;
J(n+1,n+1) = 1;
F(1) = 0;
F(n+1) = 0;
% Newton iteration
tol = h * h; % stopping tolerance
maxdiff = tol + 1; % to pass the test in first iteration and enter the loop
while maxdiff > tol
    % Assemble the Jacobian matrix DF(w) and the vector F(w)
    for i = 2 : n
        wpi = hinv2 * (w(i+1) - w(i-1));
        Pi = h2 * fz(x(i),w(i),wpi);
        J(i,i-1) = -(1 + Pi);
        J(i,i) = 2 + hh * fy(x(i),w(i),wpi);
        J(i,i+1) = -(1 - Pi);
        F(i) = 2 * w(i) - w(i-1) - w(i+1) + hh * f(x(i),w(i),wpi);
    end;
    % Solve the resulting linear algebraic equation
    s = J \ F';
    maxdiff = max(abs(s));
    w = w - s;
end;