function fig216

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]);

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

% parameters for calculation
nx = 6000;
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);
		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
grid on
box 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.02 -2 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=0.001;
g=(-y*(1-y)*z+x*y)/ep;

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

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