%
%  code27x.m -- Xtreme boundary-layer code (djm - 03 dec 01)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global alpha epsil
%alpha = 0.25;  epsil = 0.0001;
alpha = 0.25;  epsil = 0.001;

%  ODE file, max time, method & options
Xode = 'code27a';  T = 1; T0 = 0.0;  method = 'ode15s';  tol = 1e-10;

options = odeset('reltol',tol,'abstol',tol,'refine',1, ...
		'maxstep',T/20,'stats','off','events','off');
		
%  first solve
fact = (exp(1/(1-alpha))-1)/gamma(1/(1+alpha)) ...
     / (epsil^(1/(1+alpha))) * ((1+alpha)^(alpha/(1+alpha)));
IV1 = [1 fact];
[t1 y1] = feval(method,Xode,[T0 T],IV1,options);


%  second solve
IV2 = [1 0.5*fact];
[t2 y2] = feval(method,Xode,[T0 T],IV2,options);


%  final solve
y1e = y1(end,1);  y2e = y2(end,1);  ratio = (1-y2e)/(y1e-y2e);
IV3 = [1 fact*( ratio + (1-ratio)*0.5 )];
[t3 y3] = feval(method,Xode,[T0 T],IV3,options);

num2str(y3(end,1),8)

%  outer solution

y0 = exp(- (t3.^(1-alpha) - 1)/(1-alpha));

%  boundary-layer solution

yb = 1 + (exp(1/(1-alpha)) - 1) ...
   * gammainc( (t3.^(1+alpha)) / epsil /(1+alpha) ,1/(1+alpha));

figure(1);  clf
subplot(2,1,1)
plot(t1,y1(:,1),'k--');  hold on
plot(t2,y2(:,1),'k--')
plot(t3,y3(:,1),'r')

title(['\bf boundary-layer solutions (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf y(x;\epsilon=' num2str(epsil) ')'])
axis([0 0.1 0 1.1*exp(1/(1-alpha))])

subplot(2,1,2)
plot(t3,y0,'g');  hold on
plot(t3,yb,'g');
plot(t3,y3(:,1),'r')

title(['\bf outer & boundary-layer solutions (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf y(x;\epsilon=' num2str(epsil) ')'])
axis([0 0.1 0 1.1*exp(1/(1-alpha))])


figure(1);  clf
subplot(2,1,1)
plot(t1,y1(:,1),'k--');  hold on
plot(t2,y2(:,1),'k--')
plot(t3,y3(:,1),'r')

title(['\bf boundary-layer solutions (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf y(x;\epsilon=' num2str(epsil) ')'])
axis([0 1 0 1.1*exp(1/(1-alpha))])

subplot(2,1,2)
plot(t3,y0,'g');  hold on
plot(t3,yb,'g');
plot(t3,y3(:,1),'r')

title(['\bf outer & boundary-layer solutions (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf y(x;\epsilon=' num2str(epsil) ')'])
axis([0 1 0 1.1*exp(1/(1-alpha))])