function z = besselj1squared ( k )

%*****************************************************************************80
%
%% besselj1squared() computes the square of BesselJ(1, BesselZero(0,k))
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 January 2016
%
%  Author:
%
%    Original C++ version by Ignace Bogaert.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Ignace Bogaert,
%    Iteration-free computation of Gauss-Legendre quadrature nodes and weights,
%    SIAM Journal on Scientific Computing,
%    Volume 36, Number 3, 2014, pages A1008-1026.
%
%  Input:
%
%    integer K, the index of the desired zero.
%
%  Output:
%
%    real Z, the value of the square of the Bessel
%    J1 function at the K-th zero of the Bessel J0 function.
% 
  j1 = [ ...
    0.269514123941916926139021992911E+00, ...
    0.115780138582203695807812836182E+00, ...
    0.0736863511364082151406476811985E+00, ...
    0.0540375731981162820417749182758E+00, ...
    0.0426614290172430912655106063495E+00, ...
    0.0352421034909961013587473033648E+00, ...
    0.0300210701030546726750888157688E+00, ...
    0.0261473914953080885904584675399E+00, ...
    0.0231591218246913922652676382178E+00, ...
    0.0207838291222678576039808057297E+00, ...
    0.0188504506693176678161056800214E+00, ...
    0.0172461575696650082995240053542E+00, ...
    0.0158935181059235978027065594287E+00, ...
    0.0147376260964721895895742982592E+00, ...
    0.0137384651453871179182880484134E+00, ...
    0.0128661817376151328791406637228E+00, ...
    0.0120980515486267975471075438497E+00, ...
    0.0114164712244916085168627222986E+00, ...
    0.0108075927911802040115547286830E+00, ...
    0.0102603729262807628110423992790E+00, ...
    0.00976589713979105054059846736696E+00 ];

  if ( 21 < k )

    x = 1.0E+00 / ( k - 0.25E+00 );
    x2 = x * x;
    z = ...
        x *       (  0.202642367284675542887758926420E+00 ...
      + x2 * x2 * ( -0.303380429711290253026202643516E-03 ...
      + x2 *      (  0.198924364245969295201137972743E-03 ...
      + x2 *      ( -0.228969902772111653038747229723E-03 ...
      + x2 *      (  0.433710719130746277915572905025E-03 ...
      + x2 *      ( -0.123632349727175414724737657367E-02 ...
      + x2 *      (  0.496101423268883102872271417616E-02 ...
      + x2 *      ( -0.266837393702323757700998557826E-01 ...
      + x2 *      (  0.185395398206345628711318848386E+00 ) ))))))));
  else
    z = j1(k);
  end

  return
end
