Reduce ODE into 1st order

Introduction

I thought I would present a way to leverage built-in functions to do this. I've known for a long time that NDSolve sets up an ODE problem as system of first-order equations, so the basic code must be in there. Apparently, though, it's not easily found. There is the tantalizingly named Internal`ProcessEquations`FirstOrderize, which sounds perfect. It takes 4 to 6 arguments, and eventually I guessed how to set up the OP's problem.

fop = Internal`ProcessEquations`FirstOrderize[Thread[f1 == 0], {t}, 1, {x1, x2}]; 
Column[fop, Dividers -> All]

The new system is given by joining the first two components of the output:

newsys = Join @@ fop[[1 ;; 2]]

The new variables are stored in

fop[[3]]
(*  {{x1, NDSolve`x1$56$1}, {x2, NDSolve`x2$84$1}}  *)

And the relationship to the original problem is given in the rules in the last component:

fop[[4]]
(*  {NDSolve`x1$56$1 -> Derivative[1][x1], NDSolve`x2$84$1 -> Derivative[1][x2]}  *)

If you don't like the NDSolve` module variables, there is another utility that can be discovered via state@"VariableTransformation" in the NDSolve state data components. There's no documentation, AFAIK, but you can generate examples by evaluating NDSolve`StateData objects. Its form is

Internal`ProcessEquations`FirstOrderReplace[expr, indepVars, n, depVars, newVars]

(I've only ever seen n == 1 in the third position in an ODE.) For example, for the OP's system,

Internal`ProcessEquations`FirstOrderReplace[
 Thread[f1 == 0], {t}, 1, {x1, x2}, {{X, XP}, {Y,YP}}]
(*
  {ff1 + d11 X[t] + b11 XP[t] + d12 Y[t] + c11 XP[t] Y[t] + b12 YP[t] + 
     c12 X[t] YP[t] + a11 XP'[t] + a12 YP'[t] == 0,
   ff2 + d21 X[t] + b21 XP[t] + d22 Y[t] + c21 XP[t] Y[t] + b22 YP[t] + 
     c22 X[t] YP[t] + a21 XP'[t] + a22 YP'[t] == 0}
*)

Note the list of lists in the new variable names. Each instance of x1 and x1' is replaced by X and XP respectively; x1'' is replaced by XP'. Likewise for the other variable x2.

General-purpose function

Here is a function that allows renaming of the variables as the OP wishes. It's a bit tricky not to rename the left hand sides with FirstOrderReplace; I do it by inactivating Derivative temporarily.

(* With arbitrary symbol renaming *)
ClearAll[firstOrderize];
Options[firstOrderize] = {"NewSymbolGenerator" -> (Unique["y"] &)};
firstOrderize[sys_, vars_, t_, OptionsPattern[]] := 
 Module[{fop, newsym, toNewVar},
  newsym = OptionValue["NewSymbolGenerator"];
  fop = Internal`ProcessEquations`FirstOrderize[sys, {t}, 1, vars];
  If[newsym === Automatic,
   (* don't rename *)
   Flatten@ fop[[1 ;; 2]],
   (* rename *)
   toNewVar = With[{newvars = MapIndexed[newsym, fop[[3]], {2}]},
     Internal`ProcessEquations`FirstOrderReplace[#, {t}, 1, vars, newvars] &];
   Flatten@ {toNewVar[fop[[1]] /. Last[fop]],
     Activate[toNewVar[Inactivate[Evaluate@fop[[2]], Derivative]] /. 
       toNewVar[fop[[4]]]]}
   ]
  ]

Examples

OP's Example: The automatic renaming function uses Unique["y"], which will add a number to "y", whatever number is next.

firstOrderize[Thread[f1 == 0], {x1, x2}, t]
(*
  {ff1 + d11 y3[t] + b11 y4[t] + d12 y5[t] + c11 y4[t] y5[t] + 
     b12 y6[t] + c12 y3[t] y6[t] + a11 y4'[t] + a12 y6'[t] == 0, 
   ff2 + d21 y3[t] + b21 y4[t] + d22 y5[t] + c21 y4[t] y5[t] + 
     b22 y6[t] + c22 y3[t] y6[t] + a21 y4'[t] + a22 y6'[t] == 0,
   y3'[t] == y4[t], 
   y5'[t] == y6[t]}
*)

One can use the option "NewSymbolGenerator" to specify how you want the symbols generated. It should be a function, which will be applied to the NDSolve variables in fop[[3]] with MapIndexed at level {2}.

firstOrderize[{x1'[t]^2 == x2'[t] + x1[t], x2''[t] == -x1[t]}, {x1, x2}, t,
 "NewSymbolGenerator" -> (Symbol[{"a", "b"}[[First@#2]] <> ToString@Last@#2] &)]
(*  {a1'[t]^2 == a1[t] + b2[t], b2'[t] == -a1[t], b1'[t] == b2[t]}  *)

One way to get numbering from 1 to 4 every time:

Module[{n = 0},
 firstOrderize[{x1''[t] == x2[t] + x1[t], x2''[t] == -x1[t]}, {x1, x2}, t,
  "NewSymbolGenerator" -> (Symbol["y" <> ToString[++n]] &)]
 ] // Sort
(*
  {y1'[t] == y2[t],
   y2'[t] == y1[t] + y3[t],
   y3'[t] == y4[t],
   y4'[t] == -y1[t]}
*)

Here's my approach, I think it's tidier and more general:

ClearAll@to1storder

Options[to1storder] = {"form" -> (#[#2] &)};

to1storder[eq_List, func_List, argu_, OptionsPattern[]] := 
 Module[{maxorder, mapthread, lhsae, lhsde, detoae}, 
  maxorder = #[[First@Ordering[#, -1]]] &@
      Union@Cases[eq, Derivative[i_][#][_] :> i, Infinity] & /@ func;      
  mapthread = MapThread[#, {func, maxorder}] &;
  lhsae = mapthread[#[#2]@argu &];
  lhsde = mapthread[#[#2 - 1]'@argu &];
  detoae = ((f : Alternatives @@ func) (i_: 0) | Derivative[i_][f_])[a_] :> f[i][a];

  {((*Solve[*)eq /. detoae(*,#]*)/. Thread[# -> lhsde](* /. Rule -> Equal*)) &@lhsae, 
    mapthread[Table[#[n - 1]'@argu == #[n]@argu, {n, 1, #2 - 1}] &]} /. (f : 
       Alternatives @@ func)[i_] :> OptionValue["form"][f, i]]

to1storder[eq_, func_, argu_, o : OptionsPattern[]] := 
  to1storder[Flatten@{eq}, Flatten@{func}, argu, o];

If you want to put the term with highest order of derivative to the left hand side of the equation, just add the code in the anotation.

You may find the definition of detoae a little confused. To understand it, just notice that b /. (n_: 0) b -> n returns 0. For more information, check the document of Default.

The following is an example that your OrderReduce can't handle correctly:

test = With[{x = x[t], y = y[t]}, {D[x, t, t, t]^2 == -((G m x)/(x^2 + y^2)^(3/2)),
   D[y, t, t] == -((G m y)/(x^2 + y^2)^(3/2))}]

to1storder[test, {x, y}, t, "form" -> Subscript] // Flatten // TableForm

I have written the following code for this problem. The function reads as inputs the right-hand-side of the equation in implicit form and a list of it's unknown functions explicitly assigning their dependency to the time variable. Any suggestion to make this more tidy or efficient is highly welcomed.

    QDim[a_, b_] := TrueQ[Length[a] == Length[b]]
OrderReduce[f__, var__] := With[{}, If[QDim[f, var] == True,
   numbeq = Length[f];
   dimsys = 2 numbeq;
   sysvar = Array[Subscript[nvar, #][t] &, {dimsys}];
   syscom = 
    Flatten@{x, Table[D[var[[i]], {t, 1}], {i, 1, numbeq}], 
      Table[D[var[[i]], {t, 2}], {i, 1, numbeq}]};
   subvar = 
    Flatten@{Table[syscom[[i]] -> sysvar[[i]], {i, 1, dimsys}], 
      Table[D[var[[i]], {t, 2}] -> 
        D[sysvar[[numbeq + i]], {t, 1}], {i, 1, numbeq}], 
      Table[D[var[[i]], {t, 1}] -> sysvar[[numbeq + i]], {i, 1, 
        numbeq}]};
   Flatten@{f /. subvar, 
     Table[D[sysvar[[i]], {t, 1}] - sysvar[[i + numbeq]], {i, 1, 
       numbeq}]}
   , Print["Error: dimensional mismatch."]]]

If I apply this to the example above,

f1 = {a11 D[x1[t], {t, 2}] + a12 D[x2[t], {t, 2}] + 
    b11 D[x1[t], {t, 1}] + b12 D[x2[t], {t, 1}] + 
    c11 D[x1[t], {t, 1}] x2[t] + c12 x1[t] D[x2[t], {t, 1}] + 
    d11 x1[t] + d12 x2[t] + ff1, 
   a21 D[x1[t], {t, 2}] + a22 D[x2[t], {t, 2}] + 
    b21 D[x1[t], {t, 1}] + b22 D[x2[t], {t, 1}] + 
    c21 D[x1[t], {t, 1}] x2[t] + c22 x1[t] D[x2[t], {t, 1}] + 
    d21 x1[t] + d22 x2[t] + ff2};
OrderReduce[f1, {x1[t], x2[t]}] // TableForm

I correctly get:

As I said, any better alternative is highly welcomed!