I have found an answer using Solve. It's not a single function call as I'd hoped, but it could probably be made into a function. Perhaps other users will have improvements to this answer.

Consider H[L_, c] := ... as defined in the question. Define a series representation for the first-argument inverse, L[H,c]:

L[H, c] = L0[H] + 1/c^2 L2[H] + 1/c^4 L4[H]

Here only the even terms are required as can be seen in the definition of H.

Now we can setup the equation H == H[L[H, c], c] and ask Mathematica to Solve for the coefficients at each order in 1/c. For convenience, to avoid imaginary numbers in the solution, I have defined a dummy variable $\varepsilon = - 2 H$.

eqs = Map[# == 0 &, CoefficientList[(-ɛ/2) - H[L[H, c], c] + O[c, Infinity]^6, 1/c]];
Print /@ Select[eqs, ! MatchQ[#, True] &];

sol = Last@Solve[eqs, {L0[H], L2[H], L4[H]}] /. ɛ -> -2 H

Here, Last was used to take the positive solution, as discussed in the question. The final result is:

L[H, c] /. sol

It can be verified that this solution and the "by-hand" solution using ReplaceRepeated agree.

I have found an answer using Solve. It's not a single function call as I'd hoped, but it could probably be made into a function. Perhaps other users will have improvements to this answer.

Consider H[L_, c] := ... as defined in the question. Define a series representation for the first-argument inverse, L[H,c]:

L[H, c] = L0[H] + 1/c^2 L2[H] + 1/c^4 L4[H]

Here only the even terms are required as can be seen in the definition of H.

Now we can setup the equation H == H[L[H, c], c] and ask Mathematica to Solve for the coefficients at each order in 1/c. For convenience, to avoid imaginary numbers in the solution, I have defined a dummy variable $\varepsilon = - 2 H$.

eqs = Map[# == 0 &, CoefficientList[(-ɛ/2) - H[L[H, c], c] + O[c, Infinity]^6, 1/c]];
Print /@ Select[eqs, ! MatchQ[#, True] &];

sol = Last@Solve[eqs, {L0[H], L2[H], L4[H]}] /. ɛ -> -2 H

Here, Last was used to take the positive solution, as discussed in the question. The final result is:

L[H, c] /. sol

It can be verified that this solution and the t"by-hand" solution using ReplaceRepeated agree.

Edit: As pointed out by @yarchik, SolveAlways is designed precisely for setting polynomial equations to zero and solving for coefficients. While I couldn't get SolveAlways to solve for L0, L2, and L4 directly, an equivalent formulation can be achieved using an Eliminate one-liner:

Solve[Eliminate[(-ɛ/2) == H[L[H, c], c] + O[c, Infinity]^6, c], {L0[H], L2[H], L4[H]}] // Simplify // ReplaceAll[ ɛ -> -2 H]