function fig223

clf

% 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', [1925 1095 573 199]);

for ic=1:2
if ic==1
	ep=0.0008;
	x0 = 0.0;	y0 = 0.5;
	x1 = 1.0;	y1 = 2;
else
	ep=0.008;
	x0 = 0.0;	y0 = 1.5;
	x1 = 1.0;	y1 = -1.5;
end

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

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

newton_iterations=counter

if ic==1
	ya1=y;
else
	ya2=y;
end;
end


% plot computed solution
plot(x,ya1,'-','LineWidth',1,'MarkerSize',7)
hold on
plot(x,ya2,'--','LineWidth',1.1,'MarkerSize',7)
grid on
box on
axis([-0.02 1.02 -2 2])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

%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); 
hold off


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

function g=fy(x,y,z, ep)
g=(1-z+3*y^2*z)/ep;

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