Latest version, with a correction: May 7, 2000.

The routine

The routine:

\(\(makePairing[cubicHomogeneous_Function] := 
        Module[{dimension, var, toCubeCombination, recombined, 
            listOfMonomials, combinationList, temp, comb, d0, b0, 
            notFullRankQ, b, d, c, m, cm, dAnd0, a, answer}, 
          dimension = Length[Last[cubicHomogeneous]]; 
          var = Table[x[i], {i, dimension}]; 
          SetAttributes[toCubeCombination, Listable]; 
          toCubeCombination[a_ + b_] := 
            toCubeCombination[a] + toCubeCombination[b]; 
          toCubeCombination[0] := 0; 
          toCubeCombination[\((num_. )\)*\((a_)\)*\((b_)\)*\((c_)\)] /; 
              NumberQ[num] && \(! NumberQ[a]\) && \(! NumberQ[
                    b]\) && \(! NumberQ[c]\) := \((num/
                    24)\)*\((a + b + c)\)^3 + \((num/
                    24)\)*\((a - b - c)\)^3 - \((num/
                    24)\)*\((a + b - c)\)^3 - \((num/24)\)*\((a - b + c)\)^3; 
          toCubeCombination[\((num_. )\)*\((a_)\)*\((b_)\)^2] /; 
              NumberQ[num] && \(! NumberQ[a]\) && \(! NumberQ[b]\) && 
                OrderedQ[{a, b}] := \((num/6)\)*\((a + b)\)^3 + \((num/
                    6)\)*\((a - b)\)^3 - \((num/3)\)*
                a^3; 
toCubeCombination[\((num_. )\)*\((a_)\)*\((b_)\)^2] /; 
              NumberQ[num] && \(! NumberQ[a]\) && \(! NumberQ[
                    b]\) && \(! OrderedQ[{a, b}]\) := \((num/
                    6)\)*\((a + b)\)^3 - \((num/6)\)*\((b - a)\)^3 - \((num/
                    3)\)*a^3; 
          toCubeCombination[\((num_. )\)*\((a_)\)^3] /; 
              NumberQ[num] && \(! NumberQ[a]\) := num*a^3; 
          recombined = 
            toCubeCombination[
              Evaluate[Expand[cubicHomogeneous @@ var - var]]]; 
          listOfMonomials = 
            Map[If[Head[#] === Plus, List @@ #, {#}] &, recombined] // 
              Flatten; 
combinationList = 
            Union[Select[listOfMonomials, 
                  MatchQ[#1, \ \((num_. )\)*\((lin_)\)^3]\  &\ ] /. \((num_. )\)*\((lin_)\)^3 :> lin /; NumberQ[num]]; 
          temp = recombined /. \ 
              Table[combinationList[\([i]\)]^3 -> comb[i], {i, 
                  Length[combinationList]}]; 
          d0 = Table[
              Coefficient[combinationList[\([i]\)], x[j]], {i, 
                Length[combinationList]}, {j, dimension}]; 
          b0 = \(-Table[
                Coefficient[temp[\([i]\)], comb[j]], {i, dimension}, {j, 
                  Length[combinationList]}]\); 
notFullRankQ[mat_] := 
            Union[Flatten[Minors[mat, dimension]]] === {0}; 
          For[i = 0, notFullRankQ[b0], \(i++\), 
            b0 = Transpose[
                Join[Transpose[
                    b0], {\(IdentityMatrix[
                        dimension]\)[\([dimension\  - \ i]\)]}]]; 
            d0 = Join[d0, {Table[0, {dimension}]}]]; 
For[
            i = 0, notFullRankQ[d0], \(i++\), 
            b0 = Transpose[Join[Transpose[b0], {Table[0, {dimension}]}]]; 
            d0 = Join[
                d0, {\(IdentityMatrix[
                      dimension]\)[\([dimension\  - \ i]\)]}]]; 
          b = b0; 
d = d0; 
Off[
            Solve::"svars"]; 
c = \ 
            Module[{c, mat}, 
              mat = Array[
                  c, {Length[b[\([1]\)]], 
                    dimension}]; 
\(mat /. \(Solve[
                      b . mat == IdentityMatrix[dimension], 
                      Flatten[mat]]\)[\([1]\)]\) /. 
                c[i_, j_] -> 0]; 
On[
            Solve::"svars"]; 
m = \ 
            Transpose[NullSpace[b]]; 
          m = m*LCM\  @@ \ Denominator[Union[Flatten[m]]]; 
          cm = Transpose[Join[Transpose[c], \ Transpose[m]]]; 
          dAnd0\  = 
            Transpose[
              Join[Transpose[d], 
                Transpose[
                  Table[0, {i, Length[d]}, {j, 
                      Length[First[m]]}]]]]; 
a = 
            dAnd0 . Inverse[
                cm]; 
answer = \((b . c == 
                  IdentityMatrix[dimension])\) && \((Union[
                    Expand[cubicHomogeneous @@ var - var + 
                        b . \((a . c . var)\)^3]] == {0})\) && \((Union[
                    Flatten[
                      a . Transpose[NullSpace[b]]]] == {0})\) && \((Union[
                    Flatten[b . Transpose[NullSpace[a]]]] == {0})\); 
          If[answer, {a, b, c}, 
            Print["Something is wrong"]]];\)\)

Examples

Sample functions to try the routine on:

f0 = 
        Function[{x, y, z, u}, {x + (y - z + u)^3, y + (z - u)^3, 
            z + u^3, u}];

f = 
        Function[{x1, x2, x3, x4, x5}, {x1 + x2*x5^2, x2 + x1^2*x5 - x4*x5^2, 
            x3 + x2^2*x5, x4 + 2  x1*x2*x5 - x3*x5^2, x5}];

g  =  
          Function[{x1,  x2,  x3,  
              x4},  {x1  +  (x3*x1  +  x4*x2)*x4,  
              x2  -  (x3*x1  +  x4*x2)*x3,  x3  +  x4^3,  x4}];

h = 
        Function[{y1, y2, y3, y4, y5}, {y1, y2, y3, y4, y5} + 
            3 {0, 0, 
                y1^2 y2 + y1 y2^2 + 2  y1 y2 y4 - 
                  2  y1^2 y5, -y1^2 y2 - y1 y2^2 - 2  y1 y2 y3 - 
                  2  y1^2  y5, -y2^2 y3 - y2^2  y4}];

f1 = 
        Function[{x1, x2, x3, x4}, {x1, x2, x3 + x1*x2*x3 - x1^2*x4, 
            x4 + x2^2*x3 - x1*x2*x4}];

result = makePairing[f];

result[[1]] // MatrixForm