%
%  code33.m -- turning point demo I
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global alpha epsil
alpha = 4.0;  epsil = 0.0001;

%  ODE file, max time, method & options
Xode = 'code33a';  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
IV1 = [1 1/sqrt(epsil)];
[t1 y1] = feval(method,Xode,[T0 T],IV1,options);

%  second solve
IV2 = [1 0.5/sqrt(epsil)];
[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 ( ratio + (1-ratio)*0.5 )/sqrt(epsil)];
[t3 y3] = feval(method,Xode,[T0 T],IV3,options);

disp(['y(1)-1=' num2str(y3(end,1)-1,8)]);

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

plot(t3,y3(:,1),'r')
title(['\bf turning point example (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf y(x;\epsilon=' num2str(epsil) ')'])

subplot(2,1,2)

tt = 0:0.02:1;
Lg = 1 - (alpha)*tt.^(alpha-1)/2 - (1/4/epsil)*(tt.^(2*alpha));

plot(tt,Lg);  axis([0 1 -1 1])
title(['\bf turning point transformation, Q(x) (alpha=' num2str(alpha) ')'])
xlabel('\bf x-axis')
ylabel(['\bf Q(x;\epsilon=' num2str(epsil) ')'])