function [J0,Y0] = bessel_frobenius(x,P)
  
% Truncated Frobenius series for two Bessel functions
  
  a2m = 1;
  nx = length(x);
  x2m = ones(1,nx);
  J0 = a2m*x2m;
  F0 = zeros(1,nx);
  lnx = log(x);
  c2 = 2/pi;
  c1 = c2*(.577215664901532 - log(2)); % Euler's constant .577...from
                                       % Wolfram web site.
  y2 = J0.*lnx + F0; % First term in Frobenius series for y2.
  Y0 = c1*J0 + c2*y2;

  x2 = x.^2;
  Hm = 0;  % Initialize harmonic numbers with Hm=0 for m=0
  for m = 1:P
    a2m = -a2m/(2*m)^2;
    Hm = Hm + 1/m; % Hm is m'th harmonic number 1 + 1/2 + 1/3...+ 1/m
    b2m = - a2m*Hm;
    x2m = x2m.*x2;
    J0 = J0 + a2m*x2m;
    F0 = F0 + b2m*x2m;
    y2 = J0.*lnx + F0;
    Y0 = c1*J0 + c2*y2;
  end