function  [X,Y1,W2,W3,W4]=nlfd114r(Funs,I,bv,N,Nmax,TOL);
%%   function  [X,Y1,W2,W3,W4]=nlfd114r(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 with built in
%%   RICHARDSON EXTRAPOLATION TO TWO LEVELS. It calls on the algorithm
%%   nlfd114 and 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.
%%    Y1 - approximation to the solution y at X using N intervals
%%    W2 - approximation to the solution y at X extrapolating with N and 2*N 
%%    W3 - approximation to the solution y at X extrapolating with 2*N and 4*N 
%%    W4 - approximation to the solution y at X extrapolating the last 2
[X,Y1]=nlfd114(Funs,I,bv,N,Nmax,TOL);
[X2,Y2]=nlfd114(Funs,I,bv,2*N+1,Nmax,TOL);
[X4,Y4]=nlfd114(Funs,I,bv,4*N+3,Nmax,TOL);
W2=(4*Y2(1:2:2*N+3)-Y1)/3;
W3=(4*Y4(1:4:4*N+5)-Y2(1:2:2*N+3))/3;
W4=(16*W3-W2)/15;
