function [x, y] = nbvp(f,fy,fz,n,bcL,bcR)

% code to solve the nonlinear BVP
%        y'' = f(x, y, y')     for xL < x < xR  
% where 
%          y(xL) = yL ,  y(xR) = yR

% Usage (example at end of file)
%     [x, y] = nbvp(@(x,y,z)f(x,y,z),@(x,y,z)fy(x,y,z),@(x,y,z)fz(x,y,z),n,bcL,bcR)

% Output: y = solution vector at given x values

% Required arguments
%   f(x,y,z): right hand side of differential equation, with z acting as 
%             the place holder for y'
%   fy(x,y,z): derivative of f with respect to y
%   fz(x,y,z): derivative of f with respect to z
%   n: number of x points to be used between xL and xR
%   bcL = [xL, yL]: values for left BC
%   bcR = [xR, yR]: values for right BC



%%% the code begins
error = 1e-10;  % stopping error for Newton's method

% start off with a linear solution
x=linspace(bcL(1),bcR(1),n+2); 
for ix=1:n+2
	y(ix) = bcL(2) + (bcR(2)-bcL(2))*x(ix);
end

dx=x(2)-x(1);
dxx = dx*dx;
err=1;
counter=0;
while err > error
	a=zeros(1,n); c=zeros(1,n); v=zeros(1,n); u=zeros(1,n);
	for j = 2:n+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z);
    end
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:n
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
    end
	vv = v(n);
	y(n+1) = y(n+1) + vv;
	err = abs(vv);
	for jj = n:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
    end
	counter=counter+1;
    if counter>20
        error('Newton is not converging; try using a larger n')
    end
end


%%% Example: y'' = -1 + 20y + tanh(y')  for 0 < x < 1  with y(0) = 3 and y(1) = 5
%  
% bcL = [0,3]; bcR = [1,5];
% [x, y] = nbvp(@(x,y,z)-1 + 20*y + tanh(z),@(x,y,z)20,@(x,y,z)sech(z)^2,100,bcL,bcR);
%%% plot solution
% plot(x,y)

% What to do if it doesn't work (aside from checking on input errors)
%  1.  The program uses a linear function that satisfies the BCs as a start-off 
%      guess (line 22). You might try a quadratic function like  
%      y(ix) = (yL*(x(ix)-xR)^2 + yR*(x(ix)-xL)^2)/(xR-xL)^2
%  2.  If the problem has a boundary or interior layer you should consider 
%      increasing  n  (don't be shy about making n very large).
%  3.  The error value on line 16 is used by Newton's method.  It is very 
%      doubtful that you will need to change this, but certainly you can if 
%      you feel that it will help.

%  Background
%  This is a simple code that solves nonlinear 2nd order BVPs. It does not 
%  use adaptive methods, but instead uses 2nd order centered differences on 
%  a uniform grid, with Newton's method used to solve the resulting nonlinear  
%  system. The specific algorithm is described in the book Intro to Numerical  
%  Methods for Differential Equations (Holmes, 2007) on pages 59-62.  

%  It has been tested on MATLAB, version R2020b
%  Aug 20, 2020