I haven't worked out the details, but this should be doable in paraboloidal coordinates, in the same way that the geodesics on an ellipsoid were computed by Jacobi using his famous ellipsoidal coordinates (invented for that very purpose).
In the notation of the Wikipedia article, take the parameters to be $A=1$ and $B=-1$; then the coordinate surface $\mu=0$ is your paraboloid $2z=y^2-x^2$, and the geodesics on any coordinate surface should be possible to integrate (more or less) explicitly, although the details might be a bit messy. The idea is to rewrite the Hamiltonian for geodesic motion, $H=(p_x^2+p_y^2+p_z^2)/2$, in the new coordinates and apply the Hamiltion–Jacobi method (the Hamilton–Jacobi equation is separable in these coordinates).
I haven't got time to look for a good reference explaining Jacobi's work right now. I'll update the answer later if I find something.
Edit: Here's what I think is the most convenient, and least error-prone, way of setting up the equations for the geodesics. Change to paraboloidal coordinates in $R^3$; I'll call them $(u_1, u_2, u_3)$ instead of $(\lambda, \mu, \nu)$. Since this is an orthogonal coordinate system, the Euclidean metric tensor is diagonal, $ds^2=\sum_{k=1}^3 h_k^2 du_k^2$, where $h_1$, $h_2$, $h_3$ are the scale factors given in the Wikipedia article. The hyperbolic paraboloid that you are interested in is the coordinate surface $u_2=0$, which is a Riemannian manifold in itself, with coordinates $(u_1,u_3)$ and metric tensor given by $ds^2 = h_1^2 du_1^2 + h_3^2 du_3^2$ (where of course $u_2=0$ should be substituted into the expressions for the scale factors). On any Riemannian manifold, the geodesic equations are the canonical Hamiltonian equations given by the Hamiltonian function $H=\frac{1}{2} g^{ij} p_i p_j$, where $g^{ij}$ is the inverse metric tensor. In this case, we get $H(u_1,u_3,p_1,p_3) = \frac{1}{2} \left( \frac{p_1^2}{h_1(u_1,u_3)^2} + \frac{p_3^2}{h_3(u_1,u_3)^2} \right)$. So just feed this function to Mathematica, and numerically integrate the equations $$\dot{u}_1 = \partial H/\partial p_1,$$ $$\dot{u}_3 = \partial H/\partial p_3,$$ $$\dot{p}_1 = -\partial H/\partial u_1,$$ $$\dot{p}_3 = -\partial H/\partial u_3,$$ with suitable initial conditions. This will give you a geodesic emanating from a given point in a given direction. The result is in terms of paraboloidal coordinates, of course, but it is trivial to express it in terms of Cartesian coordinates (for plotting) using the formulas defining the change of variables. Finding a geodesic between two given points seems more complicated; perhaps use some "shooting" algorithm?
As I wrote above, it should be possible to integrate the equations by hand, but numerical integration seems to suffice for your purposes.
For the explicit integration of the geodesics on the ellipsoid, Jacobi's own lectures are probably as good a source as anything else (if you can read German). They are available at the Internet Archive. Elliptic coordinates are described in Lecture 26, the geodesics on the ellipsoid in Lecture 28 (p. 212).
