Is there a way to calculate the zeros of $f(z,w)= w-z^{(z^w)}$? I know the zeros of $f(z,w)=w-z^w$ have an analytic form:
$$\operatorname{zero}[z,n]=-\frac{W[-\log(z),n]}{\log(z)}$$
Is there a way to compute the zeros of
$$f(z,w)=w-z^{(z^w)}$$?
 A: Without wanting to detract too much away from Gottfried's way of solving this, I'd like to point a couple of things that may be of interest to the OP: You can actually solve many transcendental equations, by defining and inverting certain more complex functions (Note that I use a slightly different notation here: solving the equation $c^{c^z}=z$, so your $z$ correponds to my given $c$) Similar to how $z=c^z$ is solved by Lambert's $W$ function which inverts the map $zc^{-z}$, when solving the equation $zc^{-z}=1$ as $z_k=\frac{W_k(-\log(c))}{-\log(c)}$. For example, you can define the map $HW$ to be the inverse of $z\exp(\log(c)\exp(\log(c)z))$, which can solve then the equation $c^{c^z}=z$ or $f_c^{(2)}(z)=z$, where the $(2)$ denotes two-fold composition of $f_c(z)$. Such inverse maps always exist, by virtue of Lagrange's Inversion Theorem, so they can calculate solutions to similar equations fast.
For example, for a given $c\in\mathbb{C}$, then a solution of $c^{c^z}=z$ will be given as:
$$z_0=\frac{HW(-\log(c);\log(c))}{-\log(c)}$$
Using the Maple code from the back of this article, for specific $c=-2-i$ this is calculated with at least 8 digits of accuracy as: $z_0\sim 0.243918+0.1945752i$. $z_0$ along with $f_c(z_0)$ will be a 2-cycle. Check: $f_c(z_0)=1.8014674-0.9762585679$, and $f_c^{(2)}(z_0)=z_0$.
In one of the later articles it is proved that such maps ($HW$) are actually multi-valued as well and are given a recursive procedure to calculate the other branches $HW_k$, $k\in\mathbb{Z}$, without resorting to numerical methods, except to improve accuracy. So, speaking generally this seems to confirm Gottfried's result that the solutions are at least countable for each $c\in\mathbb{C}$.
Edit#1:
I am updating this to note an interesting connection: When I compare Gottfried's picture for 2-periodic points below in Update #5 with the actual image of the Julia Set for $g_c(z)=c^z$, for $c=3/2\exp(\pi i/4)$, the periodic points appear to be suspiciously close to the periodic fixed points on the outside layer of the Julia Set. I include a pic here for the Julia Set for this $c$, so you can discern yourself:

It appears that there are many many more, as Gottfried's correspond only to the outer recursive layer of the Cantor Bouquet. The main greenish circular feature is the fixed point of convergence of the iterated sequence $g_c^{(\omicron p)}(z)$. All greenish circular features in the Cantor Bouquet are also fixed points, but repelers. If you unwind the Bouquet down to smaller copies of itself, the pattern repeats around the sub-bouqets. This appears to suggest that there is a continuum of such points, or, all the fixed points are indexed by $\mathbb{Z}^{\infty}$.
I also checked the algorithm of my $HW$ functions, but unfortunately I cannot make it work to pick up more solutions. The problem seems to be related to the fact that the roots of the poly are far away from Gottfried's 2-periodic solutions and the algorithm picks up a wrong root - which then feeds to Newton and produces an overflow. I will try to optimize it a bit and see if I can make it work to pick up at least the roots that Gottfried has listed in Update#5.
A: I'll go ahead and post my code to better explain what I'm doing.  I use Mathematica.  First I define a=Log[3] as aVal=Log[3].  Then I define a twice-iterated function funB[w,n,m].  Then I use the build-in iterator NestList to iterate the function 10 times with log sheets -4 and -3 starting the iteration at $w_0=1+i$:  
aVal = Log[3];
funB[w_, n_, m_] := 1/a (Log[1/a (Log[w] + 2 n Pi I)] + 2 m Pi I);
NestList[funB[#, -4, -3] &, 1. + I, 10] // MatrixForm

$$
\left(
\begin{array}{c}
 1.\, +1. i \\
 -0.451557+4.31613 i \\
 -0.442715+4.30407 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
 -0.442731+4.3041 i \\
\end{array}
\right)
$$
And the iteration quickly settles down to w=-0.44273+43141I but that's not a 2-cycle for $3^{3^w}$
A: Spend some time reviewing Yiannis' paper on solving for the roots of the complex auxiliary equation and adapting the method to my problem $w=z^{z^w}$.  The method, if I understand it correctly, computes a Taylor series for the auxiliary equation at the origin and then uses the smallest zero in absolute value as the starting point in a Newton iteration for the root of the aux equation.  This is my code in Mathematica:
    hyperW[args_, var_, nMax_] := 
  Module[{auxEqn, g, taylorF, theZeros, smallestZero, theRoot},
   (* 
    create the auxiliary equation 
   *)
   auxEqn = Fold[Exp[#1 #2] &, Exp[z], Reverse@args];
   g[z_] := z auxEqn - var;
   (*
    create a Taylor series for the aux equation centered at zero 
   *)
   taylorF[z_] := Normal@Series[g[z], {z, 0, nMax}];
   (* solve for the zeros of the taylor series *)
   theZeros = z /. NSolve[taylorF[z] == 0, z];
   (* 
    now find smallest root of the Taylor series -- not sure why
    *)
   smallestZero = theZeros[[First@Ordering[Abs /@ theZeros, 1]]];
   (*
    solve for the root of the aux equation
    *)
   theRoot = z /. FindRoot[g[z] == 0, {z, smallestZero}];
   theRoot
   ];

I've only checked if for my expression:
theZ=3/2 Exp[Pi I/4];
auxRoot=hyperW[{-Log[theZ]},Log[theZ],10]/Log[theZ]
theZ^(theZ^auxRoot)

Out[140]= 0.630349 +0.551316 I
Out[141]= 0.630349 +0.551316 I

