Two first order non-linear differential equation: Using DSolve to obtain the solutions analytically

Amazingly, this system can be solved symbolically, at least for j2 = j1. How useful the answer may be is a different issue. If all symbols are real, then system (without initial conditions) can be rewritten as

system = {p1'[x] == j1 p1[x]^2 + 2 j1 p1[x] p2[x], 
          p2'[x] == 2 j2 p1[x] p2[x] + j2 p2[x]^2}

where p1[x] = A1[x]^2 and p2[x] = A2[x]^2. (If the symbols are complex, then p1[x] = A1[x]*Conjugate[A1[x]], p2[x] = A2[x]*Conjugate[A2[x]], and j1 and j2 replaced by their real parts.) The obvious next step, DSolve[system /. j2 -> j1, {p1, p2}, x] returns unevaluated with unexpected error messages. So, instead, eliminate p2, which is straightforward.

D[Simplify[#/p1[x]], x] & /@ First[system];
% /. Rule @@ Last[system];
Solve[First[system], p2[x]] // Flatten;
eq = Simplify[%% /. %]
(* 3 j1 j2 p1[x]^3 + 2 p1''[x] == 
   2 (j1 + j2) p1[x] p1'[x] + ((2 j1 + j2) p1'[x]^2)/(j1 p1[x]) *)

which can be solved.

s = DSolve[eq /. j2 -> j1, b1[x], x][[1, 1]]
(* p1[x] -> InverseFunction[Inactive[Integrate][4/(3 (4 j1 K[1]^2 + 
   (E^((3 C[1])/8) j1^2 K[1]^3)/(4 E^((3 C[1])/8) j1^3 K[1]^5 + 
   Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3) + 
   (4 E^((3 C[1])/8) j1^3 K[1]^5 + Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 
   16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3))), {K[1], 1, #1}] &][x + C[2]] *)

Perhaps, more convenient is

s[[2, 1]] - C[2] -> s[[2, 0, 1]][s[[1]]] - C[2];
% /. {K[1], 1, p1[x]} -> {K[1], p10, p1[x]} /. C[2] -> 0
(* x -> Inactive[Integrate][4/(3 (4 j1 K[1]^2 + 
   (E^((3 C[1])/8) j1^2 K[1]^3)/(4 E^((3 C[1])/8) j1^3 K[1]^5 + 
   Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3) + 
   (4 E^((3 C[1])/8) j1^3 K[1]^5 + Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 
   16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3))), {K[1], p10, p1[x]}] *)

p2 then can be determined by inserting the p1solution into system[[1]], and the initial condition p20 applied to obtain C[1]. If j2!= j1, DSolve returns unevaluated.

Addendum: Obtaining first integral

C[1] can be derived more simply by obtaining a first integral of eq by the standard transformation valid for autonomous second order ODEs. Doing so provides some useful insight as well.

eqv = eq /. {p1''[x] -> v'[z] v[z], p1'[x] -> v[z], p1[x] -> z}
(* 3 j1 j2 z^3 + 2 v[z] v'[z] == (j1 + j2) z v[z] + ((2 j1 + j2) v[z]^2)/(j1 z) *)

eqc1 = (DSolve[eqv, v[z], z] // First // Simplify) /. {v[z] -> p1'[x], z -> p1[x]}
(* (1/(3 j2)) 2 (j1 + j2)^2 (Log[p1[x]]/j1 + 
   (Log[(2 (j1 + j2) (j1 p1[x]^2 - p1'[x]))/(j1 p1[x]^2)] - (3 j2 
   Log[-((2 (j1 + j2) (3 j1 j2 p1[x]^2 + (-2 j1 + j2) p1'[x]))/(3 j1 j2 p1[x]^2))])
   /(2 j1 - j2))/(j1 - 2 j2)) == C[1] *)

from which C[1] is obtained by setting x == 0 and applying initial conditions. To obtain the complete solution for p1[x] requires solving eqc1 for p1'[x], which can be done for j2 = j1 and, perhaps, a few other cases. (This explains why DSolve[eq, b1[x], x] is unable to obtain solutions for arbitrary choices of j2/j1.) So, consider j2 = j1.

Simplify[eqc1 /. {p1'[x] -> j1 p10^2 + 2 j1 p10  p20, p1[x] -> p10}]
Simplify[% /. j2 -> j1]
Exp[3/8 #] & /@ %
(* -((64 j1^2 (-p10 + p20)^3)/(27 p10 p20)) == E^((3 C[1])/8) *)

This procedure also provides a simple evaluation of p2[x]. Interchange p1[x] and p2[x], and p10 and p20 in the code block just before this addendum and also in the expression for C[1] just above.