function  [X,W]=nlfd114(Funs,I,bv,N,Nmax,TOL);
%%   function  [X,W]=nlfd114(Funs,I,bv,N,Nmax,TOL);
%%   Non-Linear finite-difference algorithm to solve the boundary value problem
%%         y'' = f(x,y,y'),   a<=x<=b, y(a)=alpha, y(b)=beta
%%                   I=[a,b]; bv=[alpha,beta];
%%   Algorithm 11.4 of Burden & Faires is implemented. It calls on the 
%%   the Crout reduction algorithm for triangular systems
%% INPUTS:
%%   Funs - string (enclosed in ' ') of the name of the m-file where
%%          the corresponding system is set up to give the output
%%          values of the functions [f(x,y,y'),f_y(x,y,y'),f_y'(x,y,y')] 
%%          as  row vector [f,fy,fyprime].
%%    I  - the interval as a row vector (e.g. I=[a,b]  )
%%    bv - the boundary values as a row vector e.g. bv=[alpha,beta] 
%%    N - number of (equal) intervals is N+1 with endpoints a=x0< ... <x_N+1=b
%%    Nmax - maximum number of iterations
%%    TOL - tolerence
%% OUTPUTS:
%%    X - the N+2 points at which the approximation is made.
%%    W - approximation to the solution y at X
A=zeros(N,N);
FLAG=0;
%%Step 1
h=(I(2)-I(1))/(N+1);
  W(1)=bv(1);W(N+2)=bv(2);
W(2:N+1)=bv(1)+((bv(2)-bv(1))/(I(2)-I(1)))*h*[1:N]; %Step 2
k=1; %Step 3
while  (k<=Nmax & FLAG==0);  %Step4
    %Step 5
    x=I(1)+h;
    t=(W(3)-bv(1))/(2*h);
    [f,fy,fyprime]=feval(Funs,x,W(2),t);
    A(1,1)=2 +(h^2)*fy;
    A(1,2)=-1 +(h/2)*fyprime;
    b(1,1)=-(2*W(2)-W(3)-bv(1)+(h.^2)*f);
    %% Step 6
    for i=2:N-1;
       x=x+h;
       t=(W(i+2)-W(i))/(2*h);
       [f,fy,fyprime]=feval(Funs,x,W(i+1),t);
       A(i,i)=2 +(h^2)*fy;
       A(i,i+1)=-1 +(h/2)*fyprime;
       A(i,i-1)=-1 -(h/2)*fyprime;
       b(i,1)=-(2*W(i+1)-W(i+2)-W(i)+(h^2)*f);
    end;
    %% Step 7
       x=I(2)-h;
       t=(bv(2)-W(N))/(2*h);
       [f,fy,fyprime]=feval(Funs,x,W(N+1),t);
       A(N,N)=2 +(h^2)*fy;
       A(N,N-1)=-1 -(h/2)*fyprime;
       b(N,1)=-(2*W(N+1)-W(N)-bv(2)+(h.^2)*f);
    %% Steps 8-14
    v=crout6_7(A,b);
    W(2:N+1)=W(2:N+1)+v';
    %% Steps 15 - 17
    if norm(v,inf)<TOL;
      X=I(1)+h*[0:N+1]';W=W';
      FLAG=1;
    end;
    k=k+1;
end;
if FLAG==0;
   fprintf('Maximum number of iterations exceeded')
end;
