function [ n_data_new, n, x, fx ] = cheby_u_values ( n_data )

%% CHEBY_U_VALUES returns values of the Chebyshev U polynomial.
%
%  Differential equation:
%
%    (1-X*X) Y'' - 3 X Y' + N (N+2) Y = 0
%
%  First terms:
%
%    U(0)(X) =   1
%    U(1)(X) =   2 X
%    U(2)(X) =   4 X**2 -   1
%    U(3)(X) =   8 X**3 -   4 X
%    U(4)(X) =  16 X**4 -  12 X**2 +  1
%    U(5)(X) =  32 X**5 -  32 X**3 +  6 X
%    U(6)(X) =  64 X**6 -  80 X**4 + 24 X**2 - 1
%    U(7)(X) = 128 X**7 - 192 X**5 + 80 X**3 - 8X
%
%  Recursion:
%
%    U(0)(X) = 1,
%    U(1)(X) = 2 * X,
%    U(N)(X) = 2 * X * U(N-1)(X) - U(N-2)(X)
%
%  Norm:
%
%    Integral ( -1 <= X <= 1 ) ( 1 - X**2 ) * U(N)(X)**2 dX = PI/2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, the order of the function.
%
%    Output, real X, the point where the function is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 13;
  fx_vec = [ ...
     1.0000000000E+00,  0.4000000000E+00, -0.8400000000E+00, ...
    -0.7360000000E+00,  0.5456000000E+00,  0.9542400000E+00, ...
    -0.1639040000E+00, -1.0198016000E+00, -0.2440166400E+00, ...
     0.9221949440E+00,  0.6128946176E+00, -0.6770370970E+00, ...
    -0.8837094564E+00 ];
  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10, 11, ...
    12 ];
  x_vec = [ ...
    0.2E+00,  0.2E+00,  0.2E+00, ...
    0.2E+00,  0.2E+00,  0.2E+00, ...
    0.2E+00,  0.2E+00,  0.2E+00, ...
    0.2E+00,  0.2E+00,  0.2E+00, ...
    0.2E+00 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    x = 0.0E+00;
    fx = 0.0E+00;
  else
    n = n_vec(n_data_new);
    x = x_vec(n_data_new);
    fx = fx_vec(n_data_new);
  end

  return
end