%
%  code14.m -- IVP check (djm - 03 oct 01)
%

ep = 10.^(-2.5:0.1:1);  tt = 1.0;  tol = 1e-8;

%  evaluate exact & simple asymptotics

y_ex = log(tan(tt*ep/2 + pi/4))./ep;

y_si = tt + (ep.^2)*tt^3/6 + (ep.^4)*tt^5/24;

%  evaluate modified asymptotics

y_00 = (sqrt(2)./ep).*tan(ep*tt/sqrt(2));
y_01 = zeros(size(y_00));

for j=1:size(ep,2)
  y_01(j) = quadl('code14a',0,tt,tol,0,ep(j));
end

y_mo = y_00 + (y_01.*cosh(ep.*y_00));

%  error plot

clf;  subplot(2,1,1)

loglog(ep,abs(y_ex-y_si),'r.','markersize',10);
hold on
loglog(ep,abs(y_ex-y_mo),'b.','markersize',10);
loglog(ep,abs(y_ex-y_si),'r');
loglog(ep,abs(y_ex-y_mo),'b');  
axis([10^(-2.5) 10 1e-16 1e+5])
plot([1e-2 1],[1e-12 1e0],'k--')  
text(10^(-1.5),1e-5,'\bf slope = 6')
plot([1e-2 1],[1e-20 1e-4],'k--')  
text(10^(-0.4),1e-10,'\bf slope = 8')
text(10^(-2.1),1e3,'\bf asymptotics:  simple (red) & modified (blue)')

title('\bf asymptotic error at t=1')
  
xlabel('\bf \epsilon-axis')
ylabel('\bf |y_{ex}(1;\epsilon)-y_{as}(1;\epsilon)|')

%  uniformity plot

ep = 0.01;  dt = 0.005;  tt = dt:dt:1-dt;

%  evaluate exact & simple asymptotics

t_ex = tt*(pi/2/ep);
y_ex = log(tan(t_ex*ep/2 + pi/4))/ep;

t_si = tt*(pi/ep);
y_si = t_si + (ep^2)*(t_si.^3)/6 + (ep^4)*(t_si.^5)/24;

%  evaluate modified asymptotics

t_mo = tt*(pi/sqrt(2)/ep);
y_00 = (sqrt(2)/ep)*tan(ep*t_mo/sqrt(2));
y_01 = zeros(size(y_00));

for j = 2:size(t_mo,2)
  y_01(j) = y_01(j-1) + quadl('code14a',t_mo(j-1),t_mo(j),tol,0,ep);
end

y_mo = y_00 + cosh(ep*y_00).*y_01;

subplot(2,1,2)

plot(t_ex,y_ex,'k');
hold on
plot(t_si,y_si,'r');
plot(t_mo,y_mo,'b');
plot(pi/2/ep*[1 1],[0 1000],'k--')  
plot(pi/sqrt(2)/ep*[1 1],[0 1000],'b--')  

title('\bf non-uniform growth (\epsilon = 0.01)')
  
xlabel('\bf t-axis')
ylabel('\bf y-axis')
text(5,800,'\bf exact (black), simple (red) & modified (blue)')

axis([0 250 0 1000])