Say I have a system of ODEs $$ \dot{x} = f(x,y,t), \\ \dot{y} = g(x,y,t). $$ In this case given by the stream function $\psi(x,y,t)$ i.e $(\dot{x}, \dot{y}) = (-\partial_y\psi, \partial_x\psi) = (f, g)$, and we know the trajectories follow periodic orbits.
What numerical method would result in stable periodic orbits, if I were to plot said trajectories?
For such a system of equations, which to me looks Hamiltonian (maybe the explicit dependence of $\psi$ on $t$ could make things a bit more difficult), the best class of methods is the one of symplectic integrators, which are a subfamily of the Runge-Kutta methods. Most of them are implicit, but if your Hamiltonian $\psi$ has a specific structure (i.e., separable), then you have explicit alternatives.
In particular, symplectic methods have the favorable property of (almost) conserving the energy over long time spans, thus resulting in orbits which are stable over time. The standard explicit Euler will make your energy explose, and the implicit Euler will make it implode.
Specifically, a symplectic method of order 1 is the so-called symplectic Euler method. For order 2, you have a wider choice. For example, implicit midpoint is symplectic, as well as the Störmer-Verlet method. The best reference around for symplectic methods (and geometric integration in general) is the book Geometric numerical integration: Structure-preserving algorithms for Ordinary Differential Equations by E. Hairer, C. Lubich and G. Wanner (Springer). Chapter VI is dedicated exclusively to the problem you're dealing with.