function B = hankel_inverse ( n, x )

%*****************************************************************************80
%
%% hankel_inverse computes the inverse of a Hankel matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Miroslav Fiedler,
%    Hankel and Loewner Matrices
%
%  Input:
%
%    integer N, the order of the Hankel matrix.
%
%    real X(2*N-1), the vector that defines the matrix.
%
%  Output:
%
%    real B(N,N), the inverse matrix.
%
  x = x(:);
%
%  Define the matrix.
%
  A = hankel_matrix ( n, x );
%
%  Solve two linear systems.
%
  p = [ x(n+1:2*n-1); 0.0 ];
  u = A \ p;

  q = zeros ( n, 1 );
  q(n) = 1.0;
  v = A \ q;
%
%  Construct four matrices.
%
  z1 = zeros ( n, 1 );
  w1 = [ v(2:n); z1 ];
  M1 = hankel_matrix ( n, w1 );

  z2 = zeros ( n-1, 1 );
  w2 = [ z2; u ];
  M2 = toeplitz_matrix ( n, w2 );

  z3 = zeros ( n, 1 );
  z3(1) = -1.0;
  w3 = [ u(2:n); z3 ];
  M3 = hankel_matrix ( n, w3 );

  z4 = zeros ( n-1, 1 );
  w4 = [ z4; v  ];
  M4 = toeplitz_matrix ( n, w4 );
%
%  Construct the inverse matrix.
%
  B = M1 * M2 - M3 * M4;

  return
end

