function nbvp

% code to solve, and then plot, the solution of the nonlinear BVP
%        y'' = f(x, y, y')     for xL < x < xr  
% where 
%          y(xl) = yL ,  y(xR) = yR

% an explanation of how to use this code is given at the end of this file

% set boundary conditions
	xL=0; yL=3;
	xR=1; yR=5;

% parameters for calculation
nx = 200;
error = 0.000001;

% start off with a linear solution
x=linspace(xL,xR,nx+2); 
for ix=1:nx+2
	y(ix) = yL + (yR-yL)*x(ix);
end;


%clf
% get(gcf)
%set(gcf,'Position', [1925 1095 573 199]);

dx=x(2)-x(1);
dxx = dx*dx;
err=1;

counter=0;
while err > error

	% calculate the coefficients of finite difference equation
	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
	for j = 2:nx+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;

	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:nx
		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(nx);
	y(nx+1) = y(nx+1) + vv;
	err = abs(vv);
	for jj = nx:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;
	counter=counter+1;
end;

newton_iterations=counter

% plot computed solution
plot(x,y)
hold on
grid on
box on
%axis([-0.02 1.02 -2 2])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);
hold off

%  right hand side of differential equation
function g=f(x,y,z)
g = -1 + 20*y - tanh(z);

%  fy(x,y,z) = diff(f,y)
function g=fy(x,y,z)
g = 20;

%  fz(x,y,z) = diff(f,z)
function g=fz(x,y,z)
g = - sech(z)^2;



%  This code is run by simply entering the command  nbvp  in MATLAB (making 
%  sure that this file is in the path for MATLAB).  The user needs to first 
%  enter certain information about the problem into this file before running it.

%  Needed Input From User

%  Boundary conditions (BCs):  xL, yL, xR and yR (see lines 11-12)
%  Number of grid points:  nx (line 15)
%    
%  f(x,y,z):  this is the right hand side of the differential equation 
%             with z acting as the place holder for y' (line 80)
%  fy(x,y,z):  this is the derivative of f with respect to y (line 84)
%  fz(x,y,z):  this is the derivative of f with respect to z (line 88)

%  Example:   y'' = -1 + 20y + tanh(y')  for 0 < x< 1  with y(0) = 3 and y(1) = 5
%    nx = 200
%    xL = 0, yL = 3, xR = 1, yR = 5
%    f(x,y,z) = -1 + 20*y + tanh(z)
%    fy(x,y,z) = 20 
%    fz(x,y,z) = sech(z)^2

% 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 21). 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  nx  (don't be shy about making nx 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 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 is very robust, 
%  and was used to solve all of the nonlinear BVPs in Intro to Perturbation 
%  Methods (Holmes, 1995).  It is also relatively fast, and as an illustration 
%  is solves the above example in 0.01 sec on a laptop (and it takes only 
%  0.8 sec when nx=20000).

%  It has been tested on MATLAB, version R2009b
%  April 7, 2010