Here is an algorithm (horribly inefficient) to generate all non-hyperelliptic, non-rational, separable subfields of a non-hyperelliptic function field $F$ over a finite field $K$. Let $\Omega$ be the space of global holomorphic differentials of $F/K$. For any $K$-subspace $V$ of $\Omega$, choose a basis $v_1,\ldots,v_m$ of $V$, compute the elements $v_j/v_1,j>1$ of $F$ (and compute the algebraic relations among these $v_j$), let $E_V$ be the subfield they generate. If $E_V \ne F$ and is not rational, then you found a subfield as above. All such subfields will appear this way (proof left to the reader). There are only finitely many such $V$ since $K$ is assume finite.
Don't even dream of implementing this algorithm as is. Using the numerator of the zeta function, its factors and the Cartier operator, you can perhaps cut down the number of $V$'s that need to be tested. Maybe hyperelliptic subfields can be dealt with by using quadratic differentials.
If Florian Hess can't do it, you are probably out of luck, as far as implementation goes.
Added later: For a hyperelliptic subfield of genus $>1$, one still has a subspace $V$ but the corresponding $E_V$ is the canonical rational subfield of the hyperelliptic field. In this case, the field will be intermediate between $F$ and $E_V$ and perhaps the suggestion of Dror Speiser of using number field arguments might lead to it. It's the elliptic fields that are going to be hard to get.
