%
%  code26.m -- the next boundary layer (djm - 14 nov 01)
%

figure(1);  clf;  clear

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

epsil = 0.02;

%  evaluate exact solution:  ep*u_xx + u_x = f(x) = -exp(1-x)

u_p = (1/(1-epsil)) * (exp(1-xx)-1);
u_h = (1/(1-epsil)) * (exp(2   )-1) ...
    * (exp(-xx/epsil)-exp(-1/epsil)) / (exp(1/epsil)-exp(-1/epsil));

u_ex = u_p - u_h;

%  evaluate left BL solutions

zz = (xx+1)/epsil;

Cl0 = (exp(2)-1);  
Cl1 = Cl0 + (epsil  )*exp(2);  
Cl2 = Cl1 + (epsil^2)*exp(2);

UL0 = Cl0*(1 - exp(-zz));
UL1 = Cl1*(1 - exp(-zz)) - exp(2)*epsil*zz;
UL2 = Cl2*(1 - exp(-zz)) - exp(2)*epsil*zz ...
    + exp(2)*(epsil^2)*((zz.^2)/2 -zz);

figure(1);  clf;  subplot(2,1,1)
plot(xx,u_ex,'b');  hold on
plot(xx,u_p,'r');
plot(xx,UL0,'g:')
plot(xx,UL1,'g--')
plot(xx,UL2,'g')

title(['\bf \epsilon u_{xx} + u_x = f(x) ; (\epsilon = ' num2str(epsil) ')'])  
xlabel('\bf x-axis')
ylabel('\bf u_{exact}(x) = blue ; u_{outer}(x) = red')
axis([-1 1 0 7])

subplot(2,1,2)
plot(xx,UL0-u_ex,'g:');  hold on
plot(xx,UL1-u_ex,'g--')
plot(xx,UL2-u_ex,'g')

title(['\bf boundary-layer errors ; (\epsilon = ' num2str(epsil) ')'])  
xlabel('\bf x-axis')
ylabel('\bf u_{inner}(x) - u_{exact}(x) = green')
axis([-1 -1+10*epsil -0.1 0.1])

%  evaluate exact solution:  ep*u_xx + u_x = f(x) = -exp(1-x)

u_p = (1/(-1-epsil)) * (exp(1-xx)-exp(2));
u_h = (1/(-1-epsil)) * (1        -exp(2)) ...
    * (exp(xx/epsil)-exp(-1/epsil)) / (exp(1/epsil)-exp(-1/epsil));

u_ex = u_p - u_h;

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

title(['\bf \epsilon u_{xx} - u_x = f(x) ; (\epsilon = ' num2str(epsil) ')'])  
xlabel('\bf x-axis')
ylabel('\bf u_{exact}(x) = blue ; u_{outer}(x) = red')
%axis([-1 1 -7 0])