function fig214

% code to solve the BVP
%  y'' = f(x, y, y')  for  x0 < x < x1   '
%  y(x0) = y0  ,    y(x1) = y1

% get(gcf)
set(gcf,'Position', [654 548 573 199]);

% set boundary conditions
x0 = 0.0;	y0 = 1.0;
x1 = 1.0;	y1 = -1.0;

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

% start off with a linear solution
x=linspace(x0,x1,nx+2); 
gam=1;
for ix=1:nx+2
%	y(ix)=-gam*x(ix)*(1-x(ix));
	y(ix)=1-2*x(ix);
end;
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,'--','LineWidth',1,'MarkerSize',7)
hold on

%set(gca,'ytick',[0  0.05 0.1 0.15]);
%set(gca,'xtick',[0  0.2 0.4 0.6 0.8 1.0]);
%axis([0 1 0 0.15])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

%plot asy solution
ep=0.1;
for ix=1:nx+2
	dd=-0.75*(2*x(ix)-1)/ep;
	ya(ix)=x(ix)+1-3/(1+exp(dd));
end;
plot(x,ya,'k','LineWidth',1)

say=['\epsilon = ',num2str(ep)];
text(0.04,-1.3,say,'FontSize',14,'FontWeight','bold')

%loc='NorthWest';
loc='NorthEast';
set(gca,'FontSize',14);
legend(' Numerical',' Composite','Location',loc);
%set(findobj(gcf,'tag','legend'),'FontSize',14); 
grid on
hold off



function g=f(x,y,z)
ep=0.1;
g=(z-1)*y/ep;

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

function g=fz(x,y,z)
ep=0.1;
g=y/ep;