function fig219

clf

nx=1000;
x=linspace(0,1,nx);

% get(gcf)
set(gcf,'Position', [802 575 573 199]);

ep=0.01;

for ix=1:nx
	t1 = erf( (x(ix) - 0.5)/sqrt(2*ep) );
	t2 = exp( - (x(ix) - 0.5)^2/(2*ep) );
	ya(ix) = ( x(ix) - 0.5 )*( 1 + 5*t1 ) + 5*t2*sqrt(2*ep/pi);
	if x(ix)<0.5
		y2(ix)=-4*(x(ix)-0.5);
	else
		y2(ix)=6*(x(ix)-0.5);
	end
end;

%  y'' + p(x)y' + q(x)y= f(x)   for xL < x < xR
% set boundary conditions
	xL=0; yL=2;
	xR=1; yR=3;

h=x(2)-x(1);
% calculate the coefficients of finite difference equation
a=zeros(1,nx-2); b=zeros(1,nx-2); c=zeros(1,nx-2);
	for i=1:nx-2
		a(i)=-2+h*h*q(x(i+1),ep);
		b(i)=1-0.5*h*p(x(i+1),ep);
		c(i)=1+0.5*h*p(x(i+1),ep);
		f(i)=h*h*rhs(x(i+1),ep);
	end;
f(1)=f(1)-yL*b(1);
f(nx-2)=f(nx-2)-yR*c(nx-2);
% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);
y=[yL, y, yR];

plot(x,y,'--','Linewidth',1)
hold on
plot(x,ya,'-','Linewidth',1)
%plot(x,y2,'-.','Linewidth',1)

say=['\epsilon = ',num2str(ep)];
text(0.03,0.25,say,'FontSize',14,'FontWeight','bold')
%text(0.82,1.2,say,'FontSize',14,'FontWeight','bold')

box on
grid on
axis([0 1 -0.2 3])
loc='North';
%loc='South';

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

set(gca,'FontSize',14);
legend(' Numerical',' Composite','Location',loc);
set(findobj(gcf,'tag','legend'),'FontSize',14); 

function g=q(x,ep)
%pp=x+0.5;
pp=exp(5*(2*x-1));
g=-pp/ep;

function g=p(x,ep)
%pp=x+0.5;
pp=exp(5*(2*x-1));
g=pp*(x-0.5)/ep;

function g=rhs(x,ep)
g=0;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end