This particular problem can be solved in Maple by simple enumeration  on the set of all partitions of the original set of 8 golfers:

S:=combinat[setpartition]({$ 1..8}, 2):

S1:={S[1]}: n:=nops(S):

print(S[1]);

for i from 2 to n do

S2:=op(S1);

if `intersect`(S[i], `union`(S2))={} then

S1:={S2,S[i]}; print(S[i]) fi;

od:

Edited.

Example:

A:=Matrix(10, (i,j)->evalf(sin(i*j)^2)):

plots[matrixplot](A, axes=boxed);

In fact, you have a system of algebraic equations rather than of differential equations (no derivatives). Therefore, decide it by solve  command in real domain:

RealDomain[solve]({a^2*b=-7, a=3*b}, {a,b});

In such examples, I usually use the commands from the packages plots and plottools, as the syntax will be shorter :

restart;

sys := [x-2*y+z = 0, 2*y-8*z = 8, -4*x+5*y+9*z = -9]:

P := rhs~(op(solve(sys, [x, y, z])));  #Coordinates of the point of intersection

A := plots[implicitplot3d](sys, op([x, y, z]=~[seq(P[i]-5 .. P[i]+5, i = 1 .. 3)]), color = [blue, green, yellow]):  # Three planes

B := plottools[sphere]([op(P)], 0.3, color = red):  # The point of intersection as a small red sphere

plots[display](A, B);  # All together

There are many ways to do it. Here's one simple way:

P:=N->combinat[permute]([0$N, 1$N], N):

Example of use:

P(3);

           [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]]

One way:

restart; 

L:=[ [ [1,2], [2,1] ], [ [3,4], [4,3] ], [ [5] ] ]:

n:=0:

for i in L[1] do

for j in L[2] do

for k in L[3] do

n:=n+1: P[n]:=[op(i), op(j), op(k)]

od: od: od:

convert(P, list);

                    [[1, 2, 3, 4, 5], [1, 2, 4, 3, 5], [2, 1, 3, 4, 5], [2, 1, 4, 3, 5]]

Another more general way:

restart;

L:=[ [ [1,2], [2,1] ], [ [3,4], [4,3] ], [ [5] ] ]:

T:=combinat[cartprod](L):

n:=0:

while not T[finished]  do n:=n+1; P[n]:=op~(T[nextvalue]()) end do:

convert(P, list);                     

                       [[1, 2, 3, 4, 5], [1, 2, 4, 3, 5], [2, 1, 3, 4, 5], [2, 1, 4, 3, 5]]

Edited. 

An example:

P:=a*b^2 + b*c - a^3*c^2;

Q:=expand(P^7);

ind:=indets(%);

n:=igcd(seq(degree(Q,s), s=ind));

RealDomain[simplify](Q^(1/n), symbolic);

The example:

V:=Vector([x+y=1, x-y=2]);

A, v:=LinearAlgebra[GenerateMatrix](convert(V, list), [x,y]);  #  A is the matrix, v is the vector

<A | v>;   #  augmented matrix of the system

assume(a<>1, a*c<>1):

 is(a^2*c-a*c-a+1<>0);

                   true

Your function has an infinite number of points of inflection. To find them within a given range, we can use  RootFinding[Analytic]  command. I took the range  -5..5 .

f := x->(7-x)*sin(x^2-7):

F := diff(f(x), x, x);

L := [RootFinding[Analytic](F, re = -5 .. 5, im = -0.1 .. 0.1)];

P := map(t->[t, f(t)], L);

Visualization:

A := plot(P, style = point, symbol = solidcircle, color = red):

B := plot(f, -5 .. 5, color = blue):

plots[display](A, B);

Probably there is no need for such a large number of Digits. I left 5 digits, the expression converted to the fraction and simplified by  normal  command:

A:=2.560000000*10^(-30)*k1*(2.309486127*10^38*n-1.154743064*10^38-1.186994552*10^37*k2^2*n^2+2.373989104*10^37*k2^2*n+8.541613702*10^37*k2*n^2-1.154743064*10^38*n^2+7.119519043*10^37*n^2*k1^2+6.495962587*10^37*k1*n^2-1.293058266*10^38*k2*n^2*k1-1.186994552*10^37*k2^2-1.708322740*10^38*k2*n+1.293058266*10^38*k2*n*k1+8.541613702*10^37*k2-3.529088273*10^37*k1-2.966874313*10^37*k1*n)/(1.620847396*10^9*n-1.620847396*10^9):

normal(convert(evalf[5](A), fraction));

If you want to keep all the digits, just

normal(convert(A, fraction));

a:=Matrix([[1,2],[3,4]]):

add(add(a[i,j]*f[parse(cat(i,j))], j=1..2), i=1..2);

@vv

you are right. My approach above works only if the sets of actual  variables of the polynomials are the same. The following procedure  EquateCoeffs  solves the problem - it equates the coefficients of any two multivariate polynomials. The procedure includes the sub-procedure  coefff  that is useful in itself.

EquateCoeffs:=proc(L1::list,L2::list)

local coefff, s;

coefff:=proc(P,T,t)   # Returns the coefficient of the monomial t in the polynomial P (T - the set of polynomial's variables)

local L,H,i,k:

L:=[coeffs(P,T,'h')]: H:=[h]: k:=0:

for i from 1 to nops(H) do

if H[i]=t then k:=L[i] fi:

od:

k;

end proc:

coeffs(op(L1),'s1');

coeffs(op(L2),'s2');

s:=`union`({s1},{s2});

seq(coefff(op(L1),t)=coefff(op(L2),t), t=s);

end proc:

Example of use:

EquateCoeffs([a*x^2+b*y^2+c*x*y+f, [x,y]], [d*x^2+e, [x]]);

                              f = e, a = d, b = 0, c = 0

restart;

P1:=(1/2*(p*a*b+(g-p)*b-g))*b*v*a*ln(E)^2-(-1+b)*v*(g-p+a*p)*b*a*ln(E)*ln(K)-b*p*(a-1)*v*a*ln(E)*ln(L)+v*a*b*ln(E)+(1/2*(p*(-1+b)*a+(g-p)*b+p))*(-1+b)*v*a*ln(K)^2+(-1+b)*v*p*(a-1)*a*ln(K)*ln(L)-v*a*(-1+b)*ln(K)+(1/2)*a*p*v*(a-1)*ln(L)^2-v*(a-1)*ln(L):

P2:=x_1*ln(E)+x_11*ln(E)^2+x_12*ln(E)*ln(K)+x_13*ln(E)*ln(L)+x_2*ln(K)+x_22*ln(K)^2+x_23*ln(K)*ln(L)+x_3*ln(L)+x_33*ln(L)^2:

coeffs(P1, [ln(E),ln(K),ln(L)])=~coeffs(P2, [ln(E),ln(K),ln(L)]);

In your example, it is more convenient to work with sets than with lists:

A:={x, y, x^2*y, x*y^2, y^2}:   B:={x^2, y^3}:

C:=A:

for a in A do

if convert([seq(degree(a/b, x)>=0 and degree(a/b, y)>=0, b=B)],`or`) then C:=C minus {a} fi;

od:

C;

                                        {x, y, y^2, x*y^2}

Of course, it's easy to rewrite this code as a procedure with any number of variables.