% Script to plot comparison of  Bessel functions J0(x) and K0(x) to their 
% large-x asymptotic approximations

  x = 0:0.1:10;
  nu = 0;
  A1 = sqrt(2./(pi*x)); % First-order amplitude
  phi0 = x - nu*pi/2 - pi/4; % Zeroth order phase
  phi2 =  (nu^2 - 1/4)./(2*x); % Second order correction to phase
  Jnu1 = A1 .* cos(phi0);
  Ynu1 = A1 .* sin(phi0);
  Jnu2 = A1 .* cos(phi0 + phi2);
  Ynu2 = A1 .* sin(phi0 + phi2);
  
  plot(x,Jnu1,'r-',x,Jnu2,'g-',x,besselj(0,x),'b-')
  xlabel('x')
  ylim([-2 2])
  title(['Asymptotic series approximations exp(S0+...SP) to J0(x) (solid) and' ...
	 ' Y0(x) (dashed)'])
  legend('P=1','P=2','Exact')
  hold on
  plot(x,Ynu1,'r--',x,Ynu2,'g--',x,bessely(0,x),'b--')
  hold off