function fig45b

clf

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

nx=500;

aa=1; ba= 1; ep=1/25;
t0 = sqrt(3) - log(2 + sqrt(3))/2;
a0 = cos(t0/ep - pi/4);
a1 = cos(t0/ep + pi/4);
aR = (aa*3^(1/4) - ba*a1*exp(-pi/(4*ep))) / (2*a0 - a1*exp(-pi/(2*ep)));
bR = exp(-pi/(4*ep))*(ba - aR*exp(-pi/(4*ep)));

A = 2^(5/6)*sqrt(pi)*aR;
B = bR*sqrt(pi)/2^(1/6);


xxL=linspace(-1,0.001,nx);
for ix=1:nx
	t = ( 1 - xxL(ix) )*sqrt( xxL(ix)^2 - 2*xxL(ix) ) /2 - log( 1 - xxL(ix) + sqrt( xxL(ix)^2 - 2*xxL(ix) ) )/2;

	qa = xxL(ix)*(2-xxL(ix));
	yyL(ix) = (2*aR*cos(t/ep - pi/4) + bR*cos(t/ep + pi/4))/abs(qa)^(1/4);

end;

xxR=linspace(0.001,1,nx);
for ix=1:nx
	k = (xxR(ix)-1)*sqrt(2*xxR(ix)-xxR(ix)^2)/2 + asin(xxR(ix)-1)/2 + pi/4;
	qa = xxR(ix)*(2-xxR(ix));
	yyR(ix) = (aR*exp(-k/ep) + bR*exp(k/ep))/qa^(1/4);

end;


%  y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr
% set boundary conditions
	xl=-1; yl=1;
	xr=1; yr=1;

x=linspace(-1,1,nx);
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];


hold on
box on
grid on

plot(xxl,yyl,'--b','linewidth',1.3)

plot(xxr,yyr,'-.k','linewidth',1.3)

plot(x,y,'-r','linewidth',1.3)

axis([-1 1 -3 6])
%loc='NorthWest';
loc='NorthEast';

xlabel('x-axis','fontsize',14,'fontweight','bold')
ylabel('solution','fontsize',14,'fontweight','bold')

set(gca,'fontsize',14);
legend(' y_l',' y_r',' y','location',[0.72788107277636 0.62291542458053 0.114672364672365 0.146186440677966]);
set(findobj(gcf,'tag','legend'),'fontsize',14); 


function g=q(x,ep)
g=-x*(2-x)/ep^2;

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

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