function A = sylvester_matrix ( p1, p2 )

%*****************************************************************************80
%
%% sylvester_matrix() returns the SYLVESTER matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 January 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jacqueline Burm, Paul Fishback,
%    Period-3 Orbits Via Sylvester's Theorem and Resultants,
%    Mathematics Magazine,
%    Volume 74, Number 1, February 2001, pages 47-51.
%
%  Input:
%
%    real p1(n1), p2(n2): two polynomials.
%    On input, the constant coefficient is listed first.
%
%  Output:
%
%    real A(n1+n2,n1_n2): the Sylvester matrix.
%
  n1 = length ( p1 ) - 1;
  n2 = length ( p2 ) - 1;

  p1 = p1(end:-1:1);
  p2 = p2(end:-1:1);

  A = zeros ( n1+n2, n1+n2 );

  for i = 1 : n2
    A(i,i:n1+i) = p1(1:n1+1);
  end

  for i = 1 : n1
    A(i+n2,i:n2+i) = p2(1:n2+1);
  end

  return
end
