%
%  code25.m -- my first boundary layer (djm - 14 nov 01)
%

figure(1);  clf;  clear

N=400;  xx = -1:2/N:1;

alleps = [0.01 0.02 0.04];

for j=1:3
  epsil = alleps(j);

%  evaluate exact solution

  u_1 = (1/(1-epsil^2)) * (exp(1-xx));
  u_2 = (1/(1-epsil^2)) * (- exp(2) * sinh((1-xx)/epsil)/sinh(2/epsil));
  u_3 = (1/(1-epsil^2)) * (- sinh((1+xx)/epsil)/sinh(2/epsil));

  u_ex = u_1 + u_2 + u_3;

  figure(1);  subplot(2,1,1)
  plot(xx,u_ex,'b');  hold on
  plot(xx,u_1,'r');

  title(['\bf boundary-layer examples (\epsilon = ' num2str(alleps) ')'])  
  xlabel('\bf x-axis')
  ylabel('\bf u_{exact}(x) = blue ; u_{outer}(x) = red')
  axis([-1 1 0 8])

  subplot(2,2,3)
  plot(xx,u_1-u_ex,'b');  hold on

  title(['\bf error (\epsilon = ' num2str(alleps) ')'])  
  xlabel('\bf x-axis')
  ylabel('\bf u_{outer}(x) - u_{exact}(x) = blue')
  axis([-1 -0.8 0 exp(2)])

  subplot(2,2,4)
  plot(xx,u_1-u_ex,'b');  hold on

  title(['\bf error (\epsilon = ' num2str(alleps) ')'])  
  xlabel('\bf x-axis')
  ylabel('\bf u_{outer}(x) - u_{exact}(x) = blue')
  axis([0.8 1 0 1.1])
end

disp('time for a pause')
pause

eu_xx = (epsil^2) * u_1 + u_2 + u_3;

figure(2);  clf
plot(xx,abs(u_ex),'b');  hold on
plot(xx,abs(u_1),'r--');
plot(xx,abs(eu_xx),'k');

title(['\bf termwise comparison (\epsilon = ' num2str(epsil) ')'])  
xlabel('\bf x-axis')
ylabel('\bf u = blue ; \epsilon^2 u_{xx} = black ; forcing = red')

disp('time for a pause')
pause

zz = (1+xx)/epsil;
fL = exp(2);  fpL = -fL;  fppL = fL;

u_L = fL*(1 - exp(-zz)) + fpL*epsil*zz + fppL*(epsil^2)*(2 + zz.^2)/2;
u_L = u_L + 0*0.1*exp(zz);

zz = (-1+xx)/epsil;
fR = 1     ;  fpR = -fR;  fppR = fR;

u_R = fR*(1 - exp( zz)) + fpR*epsil*zz + fppR*(epsil^2)*(2 + zz.^2)/2;

figure(3);  clf
plot(xx,u_ex,'b');  hold on
plot(xx,u_L,'g--');
plot(xx,u_R,'g--');
plot(xx,u_1,'r--');

title(['\bf boundary-layer analysis (\epsilon = ' num2str(epsil) ')'])  
xlabel('\bf x-axis')
ylabel('\bf u_{exact} = blue ; u_{outer} = red ; u_{bl} = green')
axis([-1 1 0 8])