DSolve leaks internal error messages when attempting to solve two first order autonomous ODEs

Question

Issue reported to Wolfram, Inc as a possible bug in Version 12.1.1; CASE:4630268.

With Mathematica "12.1.1 for Microsoft Windows (64-bit) (June 19, 2020)", DSolve produces unexpected error messages:

 DSolve[{p1'[x] == p1[x]^2 + 2 p1[x] p2[x], 
         p2'[x] == 2 p1[x] p2[x] + p2[x]^2}, {p1, p2}, x]

Union::normal: Nonatomic expression expected at position 2 in { ... }⋃$Failed.

Flatten::normal: Nonatomic expression expected at position 1 in Flatten[$Failed].

and returns unevaluated after a few minutes. Evidently, DSolve has passed the bad argument {...}⋃$Failed to Union. I am asking

• Is this a bug? (I think it is, because DSolve should not show internal error messages to the user.)
• Is there a work-around? (This system is solvable, when converted to a single second-order ODE, so the pair of first-order ODEs also should be solvable.)

Update: Issue resolved in Mathematica "12.3.1 for Microsoft Windows (64-bit) (June 19, 2021)"

This bug has been eliminated in the sense that the internal errors no longer are leaked. Instead, DSolve simply returns unevaluated after a few minutes.

Answer 1

Does this count as a workaround?:

(* dividing 2nd ODE by 1st yields a homogeneous ODE *)
p2sol = DSolve[
     {p2'[p1] == (2 p1 p2[p1] + p2[p1]^2) /
                  (p1^2 + 2 p1 p2[p1])}, p2, p1] /.
    C[1] -> Log[C[1]] /. 
   p_Power :> RuleCondition[p, True];

(* p2sol turns the p1'[x] ODE in the system into a separable equation *)
PrintTemporary@Dynamic[foo = Clock[Infinity]];
TimeConstrained[
   (Print[foo]; #) &@
    Flatten@ DSolve[#, p2, x],
   30,
   Print[Style[foo, Red]]; $Failed] & /@
 (Last[system] /. p1[x] -> p1[p2[x]] /. p2sol)

(*
3.30611
33.329
63.273
{{p1 -> Function[{x}, 
    InverseFunction[
      Inactive[
         Integrate][(-9 C[1] K[1]^2 + 
          Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(1/3)/(
        K[1] (-2 2^(1/3) 3^(2/3) C[1] K[1] + 
           9 K[1] (-9 C[1] K[1]^2 + 
              Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(1/3) + 
           2^(2/3) 3^(
            1/3) (-9 C[1] K[1]^2 + 
              Sqrt[3] Sqrt[C[1]^2 K[1]^3 (4 C[1] + 27 K[1])])^(
            2/3))), {K[1], 1, #1}] &][
     x/3 + C[2]]]}, $Failed, $Failed}
*)

One can combine with with p2sol to obtain p2'[x]. The solutions can be stated as implicit equations, but Mathematica tries really hard to solve them.

Note this system and the one it came from admit two one-parameter families of symmetries, scalings {p1, p2, 1/x} -> C[1] {p1, p2, 1/x} and translations x -> x + C[2]. Thus they can theoretically be expressed as successive quadratures as above, provided one can solve the intermediate equations such as the one produced by this generalization:

DSolve[
  {p2'[p1] == (2 a p1 p2[p1] + b p2[p1]^2)/(c p1^2 + 2 d p1 p2[p1])}, 
  p2, p1]

(*
Solve[(-c (b - 2 d) Log[p2[p1]/p1] +
   (b c - 4 a d) Log[-2 a + c - (b p2[p1])/p1 + (2 d p2[p1])/ p1]) / 
  ((2 a - c) (b - 2 d)) == C[1] - Log[p1], p2[p1]]
*)