function sbvp

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

%  the goal is to calculate the solution when ep = 10^eMin for eMin negative

%  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;

% number of grid points
nx = 1000;

%  number of ep values used to reach ep = 10^eMin starting from ep = 10^eMax
nep = 160;
eMax = 0;
eMin = -4;

%  error used in Newton's method
error = 0.00001;

% 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;


dx=x(2)-x(1);
dxx = dx*dx;
eps=linspace(eMin,eMax,nep);
for ip=1:nep
ep=10^eps(nep-ip+1)
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,ep);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z,ep);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z,ep);
	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;

	if counter == 100
		fprintf('\n The program has paused because there have been 100 Newton iterations for ep = %8.4e\n Its recommended that you stop execution (ctrl+c or ctrl+Break) \n and increase either neap or nx, or increase eMin \n',ep); 
		pause
	end


end;

newton_iterations=counter;

end;

% plot computed solution
plot(x,y)
hold on
grid on
box on

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
%axis([-0.05 1.05 1 5])
set(gca,'FontSize',14);
hold off


function g=f(x,y,z,ep)
g=(-tanh(z)+y-1)/ep;

function g=fy(x,y,z,ep)
g=1/ep;

function g=fz(x,y,z,ep)
g=-(sech(z))^2/ep;



%  This code is run by simply entering the command  sbvp  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 14-15)
%  Number of grid points:  nx (line 18)
%    
%  f(x,y,z,ep):  this is the right hand side of the differential equation 
%                with z acting as the place holder for y' (line 99)
%  fy(x,y,z,ep):  this is the derivative of f with respect to y (line 102)
%  fz(x,y,z,ep):  this is the derivative of f with respect to z (line 105)

%  The values of ep used in calculation are printed out.  To suppress this, 
%  simply add a  ;  to the end of line 39.

%  Example:   ep*y'' = -1 + y - tanh(y')  for 0 < x < 1  
%             with  y(0) = 3 and y(1) = 5
%    xL = 0, yL = 3, xR = 1, yR = 5
%    nx = 1000
%    nep = 160
%    eMax = 0
%    eMin = -4
%    f(x,y,z,ep) = (-1 + y - tanh(z))/ep
%    fy(x,y,z,ep) = 1/ep
%    fz(x,y,z,ep) = - sech(z)^2/ep

%  Explanation of Example: The objective is to compute the solution for 
%  ep = 10^{-4}. The code does this by starting with ep = 10^eMax = 1 and
%  then finding the solution by successively reducing ep down to 
%  ep = 10^eMin = 10^{-4} (taking nep = 160 steps to do this).  The code
%  is reasonably fast, taking about 1.2 sec on a desktop computer.

% What to do if it doesn't work (aside from checking on input errors)
%  1.  The code pauses if the number of Newton iterations, for a given value 
%      of ep, reaches 100.  If this happens, it is recommended that you stop 
%      execution and increase the value of nep (line 21), and possibly the 
%      value of nx (line 18).  For example, if eMin = -6 then the above code 
%      will pause.  Stopping and taking  nx = 10000  and  nep = 2500 
%      then all is well, and the solution is computed. 
%  2.  Given the nonlinear nature of the problem, the first initial guess 
%      for the solution, which is a linear function, is not adequate (line 31). 
%      In such a case, you might try a quadratic function like  
%      y(ix) = (yL*(x(ix)-xR)^2 + yR*(x(ix)-xL)^2)/(xR-xL)^2
%  3.  The error value on line 26 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 which have 
%  boundary or interior layers.  The reason these are difficult to solve is 
%  that the stiffness of the problem means that the nonlinear solver (e.g., 
%  Newton's method) needs a very good guess for the solution for it to be 
%  able to converge.  The smaller the value of ep, the better the guess must be.  
%  This code does this by first solving the problem for a nonstiff value of 
%  ep (e.g., ep = 1), and then using the computed solution as the starting 
%  guess for a smaller value of ep.  By successively reducing ep it becomes
%  more likely that it is possible to reach the value of ep^eMin (so this 
%  is a homotopy method for singularly perturbed BVPs).  The more negative 
%  eMin the larger you will need to take nep (line 21).  Also, you should 
%  make sure there are several grid points within the layer, so if the layer 
%  is very narrow you will need to use fairly large values for nx (line 18).  
%  The BVP solver used here is the same one as used in  nbvp.m  

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