%
%  code30.m -- gumdrop boundary-layer (28 nov 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

disp(' ')
Up = input('initial derivative (U''):  ');

%  ODE file, IVs, final time, method & options

Xode = 'code30a';  IVs = [0 0 0];  
method = 'ode45';  tol = 1e-12;
T = 2*Up + (1+log(4))/2;  

options = odeset('reltol',tol,'abstol',tol,'refine',1, ...
		'maxstep',T/20,'stats','off','events','off');

%  ODE solve with period & amplitude determination

[t y] = feval(method,Xode,[0 T],IVs,options);
Tend = t(end);

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  solution plot

figure(3);  clf;
subplot(2,1,1)

plot(t,y(:,1),'b');  hold on
yp = -2*tanh(t).*y(:,1)-log(cosh(t));
plot(t,yp,'r')
plot(t,-(t/2)+log(exp(1/4)*sqrt(2)),'g--')

title('\bf gumdrop boundary correction (with linear asymptote)')
xlabel('\bf T-axis');  ylabel('\bf Y_1(T) in blue, Y_1\prime(T) in red');

amax = max(y(:,1));    amin = min(y(:,1));  
Tend = max(t);  
axis([0 4 -2 1])

disp(' ')
disp(['slope                = ' num2str(yp(end),8)]);
disp(['intercept            = ' num2str(y(end,1)-yp(end)*Tend,8)]);
disp(['ln(exp(1/4) sqrt(2)) = ' num2str(0.25+log(2)/2,8)]);
disp(['magic constant       = ' num2str(y(end,3),8)]);
disp(' ')

subplot(2,1,2)
Yfc  = tanh(t) + (1/Up)*y(:,1);
Yfcp = (sech(t)).^2 + (1/Up)*yp;
plot(t,Yfc ,'b');  hold on
plot(t,Yfcp,'r');

title(['\bf gumdrop boundary solution to 1^{st} correction (A^2 = ' ...
      num2str(Up,8) ')'])
xlabel('\bf T-axis');  ylabel('\bf Y(T) in blue, Y\prime(T) in red');
axis([0 T -0.1 1.1])

figure(1);
subplot(1,2,2)
plot(sqrt(Up)*tanh(t),Up*sech(t).^2,'g');  hold on
plot(sqrt(Up)*Yfc,Up*Yfcp,'r')
hold off

title('\bf boundary-layer phase plane')
xlabel('\bf u(t)');  ylabel('\bf u\prime(t)');
axis equal
axis([-1.1*sqrt(Up) 1.1*sqrt(Up) -0.5 1.1*Up])