How to get series coefficients from functional equations?

Question

$p(z),t(z)$ are two mutually defined functional equations, while $\widehat{G}(z)$ is the exponential generation function of A182173 (maybe, I am not sure...lol):

$$\begin{cases} p(z)=e^{t(z)}-t(z)+2 z-1\\ t(z)=2\left(e^{p(z)}-e^{p(z)/2}+z\right)-p(z)\\ \widehat{G}(z)=p(z)+t(z)-2 z \end{cases}$$

Is there any way to get series coefficients effectively?

p[z]->2z+E^t[z]-1-t[z]
t[z]->2*( z+ E^p[z] -E^(p[z]/2) ) - p[z]
Ge=p[z]+t[z]-2z
A182173[n_]:=n!SeriesCoefficient[Ge,{z,0,n}]

Answer 1

Find Ge[z]:

r1 = p[z] -> 2 z + E^t[z] - 1 - t[z]

r2 = t[z] -> 2*(z + E^p[z] - E^(p[z]/2)) - p[z]

Ge[z_] = p[z] + t[z] - 2 z /. r1 /. r2 // FullSimplify

(*   -1 + E^(2 (-E^((p[z]/2)) + E^p[z] + z) - p[z])   *)

SeriesCoefficient produces derivatives of p at z == 0. Therefore we have to solve for the p[0], p'[0], p''[0] ...

der1[i_] := Derivative[i][p][0] == Derivative[i][2 # + E^t[#] - 1 - t[#] &][0]

der2[i_] := 
  Derivative[i][t][0] == Derivative[i][2*(# + E^p[#] - E^(p[#]/2)) - p[#] &][0]

sol[0] = First @
   Solve[{der1[0], der2[0]}, {Derivative[0][p][0], 
           Derivative[0][t][0]}, Reals]

(*   {p[0] -> 0, t[0] -> 0}   *)

Lower derivatives are needed to solve for higher ones:

sol[j_ /; j > 0] := sol[j]=
   First@Solve[{der1[j], der2[j]} //. 
   Flatten[Table[sol[k], {k, 0, j - 1}]], {Derivative[j][p][0], 
   Derivative[j][t][0]}, Reals] 

Inserting these found derivatives of p into the series coefficients gives the desired A numbers.

A[n_] := n! SeriesCoefficient[Ge[z], {z, 0, n}] //. 
  Flatten[Table[sol[i], {i, 0, n}]]

Table[A[i], {i, 1, 10}]

(*   {2, 10, 94, 1466, 31814, 887650, 30259198, 1218864842, 56644903958,
      2983300619410}   *)

Answer 2

As Akku14 notes in his answer, one can determine from the conditions that both functions $p,t$ have a zero constant term. With that additional piece of information, we can use SolveAlways[] to determine the rest of the required coefficients:

n = 10; (* desired expansion order *)
p = Sum[pc[k] z^k, {k, 0, n}] + O[z]^(n + 1);
t = Sum[tc[k] z^k, {k, 0, n}] + O[z]^(n + 1);

vals = SolveAlways[{p == Exp[t] - t + 2 z - 1, t == 2 (Exp[p] - Exp[p/2] + z) - p} /.
                   {pc[0] -> 0, tc[0] -> 0}, z];

Rest[CoefficientList[Sum[(pc[k] + tc[k]) z^k, {k, 0, n}] - 2 z + O[z]^(n + 1) /.
     Join[{pc[0] -> 0, tc[0] -> 0}, First[vals]], z]] Range[n]!
   {2, 10, 94, 1466, 31814, 887650, 30259198, 1218864842, 56644903958,
    2983300619410}