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

function [w x] = fpifd(f,a,b,ya,yb,n)
h = (b - a) / n;
rho = h * h * 10; % parameter of the fixed point iteration
c1 = 0.5 / (1 + rho);
c2 = rho / (1 + rho);
c3 = -h * h * 0.5 / (1 + rho);
c4 = 0.5 / h;
% The mesh points and initial guess
x = linspace(a,b,n+1);
w = linspace(ya,yb,n+1);
% Fixed point iteration
tol = h * h; % stopping tolerance
maxdiff = tol + 1; % to pass the test in first iteration and enter the loop
while maxdiff > tol
    maxdiff = 0;
    new1 = w(1);
    for i = 2 : n
        new = c1 * (w(i-1) + w(i+1)) + c2 * w(i) + c3 * f(x(i), w(i), (w(i+1) - w(i-1)) * c4);
        diff = abs(w(i) - new);
        if diff > maxdiff 
            maxdiff = diff;
        end;
        w(i-1) = new1;
        new1 = new;
    end;
    w(n) = new1;
end;