computing monodromy of branched cover of sphere specified by a polynomial We know by the Riemann Existence Theorem that any Riemann surface can arise holmorphically as the branched cover of a sphere:


*

*Which Riemann surfaces arise from the Riemann existence theorem?
Therefore, any branched cover of the sphere can be represented by a rational function polynomial $\frac{p (z)}{q(z)}: X \to \mathbb{C}$.
Given arbitrary polynomial $p$, do there exist ways to compute the monodromy around each critical points?  
For example, $z \mapsto z^2 + \frac{1}{z^2}$ (not a polynomial) is 4-sheeted cover of the Riemann sphere $\hat{\mathbb{C}}$. It has 3 branch points, where
$$p'(z) = 2z - \frac{2}{z^3} = 0 $$ 
It seems to have 4 branch points at $z = 1, i, -1, -i$.
How do we calculate of which sheets were permuted given $p(z)$?

This is motivated by trying to understand a proof of the Riemann Existence Theorem in the case of genus 0 branched covers of the sphere $S^2$.

This question is the reverse direction of my previous question on branched covers:
polynomial branched cover of the sphere with specified monodromy
Our example has a symmetry - and I am trying to avoid techniques which exploit them too much.
 A: In your example, everything is easy to calculate. Use the fact that your $f$  is a composition of $f(z)=z+1/z$ and $h(z)=z^2$. The plane is partitioned into $8$ regions
by the coordinate axes and the unit circle. Each region corresponds to a sheet
over upper or lower half-planes. The sheets are connected in the same way as the
regions in the $z$-plane.
You found the critical points incorrectly: $0$ and $\infty$ are missing.
In general, a rational function of degree $4$ has $6$ critical points.
The same method works in general (but in practice becomes complicated for functions of
higher degree). You find critical points and critical values. Then make convenient cuts
between critical values to break the sphere into simply connected regions. Then find out
how components of the preimage of these regions tile the sphere of the independent variable.
In your example, critical values are $\pm1$, $0$ and $\infty$. So it is reasonable to
cut the image sphere into upper and lower half-planes. Finding the components of preimage
of these half-planes is simplified by the fact that your function is a composition of
two simpler functions. In general, an approximate calculation is sufficient to see how the
components of preimage fit together.
EDIT. In the classical literature, they cut the image into two hemispheres by drawing
a Jordan curve passing through all critical values. This defines a cell decomposition of the image sphere. Then they take the full preimage of this cell decomposition, or the full preimage of the dual call decomposition. In the second case one obtains an embedded
bipartite graph (1-skeleton of the preimage cell decomposition). This graph is called
the line complex; it permits to visualize the monodromy. It shows you how the half-sheets
are pasted together. Of course everything depends on
the choice of the Jordan curve, and there is no canonical way to choose it.
Once it is chosen, the line complex can be computed if you can solve the equation $f(z)=a$
numerically. The line complex is what replaces the dessin d'enfant in the case of more than 3 critical values.  
A: David linked to some computational software implementing it; I just want to mention that this can also be done in Maple.  I'd rather I was referencing either Sage because it's free and open source, or Bertini because it's free and Dan Bates here at CSU is one of the main authors, but I'm not actually much of a computing guy and this is what I figured out how to do in a few minutes with  google.  
For your example, "monodromy(y*z^2-z^4-1,y,z)" returns

[-3.60000000000, [-1.81567148865360 I, -0.550760424586247 I,  0.550760424586247 I, 1.81567148865360 I], [[-2., [[1, 2], [3, 4]]],    [2., [[1, 3], [2, 4]]], [infinity, [[1, 4], [2, 3]]]]]

Where the output means the following: -3.6 is a regular value of your map, with four pure imaginary roots given.  There are three critical values, -2,2, and infinity, and taking a path from -3.6 around your the given critical point (I'm not sure exactly which path) and back to -3.6, and lifting it to the cover, permutes the sheets in the following way.
