Warning: long post ahead. To simply watch an animation, scroll to the bottom ;)
Requirements on the wavepacket
First let's do some estimations, using the exact analytic expression for evolution of a Gaussian wave packet in free space. Hopefully, it'll not be too wrong in the case when Coulomb potential is on. The expression looks like this:
$$\operatorname{gaussian}(x,t)=\sqrt{\frac{\sigma_0}{\sqrt{2\pi}}\bigg/\left(\sigma_0^2+\frac{i\hbar t}{2m}\right)}
\exp\left(-\sigma_0^2\left(\frac p\hbar\right)^2\right)
\exp\left(-\frac{\left(x-2i\sigma_0^2\frac p\hbar\right)^2}{4\left(\sigma_0^2+\frac{i\hbar t}{2m}\right)}\right).
\tag1$$
From it we can find that expected value of $x$ changes in time as
$$\langle x\rangle=\int\limits_{-\infty}^\infty \left|\operatorname{gaussian}(x,t)\right|^2 x\,\mathrm{d}t=\frac{pt}m,
\tag2$$
while standard deviation $\sigma$ of position evolves as
$$\sigma(t)=\sqrt{\langle x^2\rangle-\langle x\rangle^2}=\sqrt{\sigma_0^2+\frac{t^2\hbar^2}{4m^2\sigma_0^2}}.
\tag3$$
We'll require that after two orbital periods, $\sigma$ becomes no more than 2 times larger than initially, i.e.
$$\sigma(2T)\le2\sigma_0.
\tag4$$
If we consider motion of a classical particle in Coulomb field
$$V(r)=-\frac\alpha r,
\tag5$$
where $\alpha=\frac{-Qq}{4\pi\varepsilon_0}$ and $Qq$ is product of charges of the attracting center and the moving particle, then orbital period is$^\dagger$
$$T=2\pi a^{3/2}\sqrt{\frac m\alpha},
\tag6$$
where $a$ is semi-major axis of the elliptic orbit.
Thus, solving inequality $(4)$, we get:
$$\sigma_0\ge\sqrt{\frac{2\pi\hbar}{\sqrt{3\alpha m}}}a^{3/4}.
\tag7$$
Thus we have the lower bound on $\sigma_0$ to limit the rate of spreading of our wave packet. But we also want the wave packet to be localized enough, so that it would be easy to visually separate it from the nucleus. We'll require that
$$\sigma_0\le a/20,
\tag8$$
and then, taking $\sigma_0$ to be the lowest possible according to $(7)$, we have our lower bound on $a$:
$$a\ge\frac{(800\pi\hbar)^2}{3\alpha m}.
\tag9$$
Explicit example of working parameters
OK, let's go numeric now. Based on the above estimation, we'll choose some parameters for electron's orbit and "experimentally" check whether the estimations yield anything good.
In general, for a hydrogen-like system we have $Qq=-Ze^2$, where $Z$ is the atomic number of the nucleus. Inserting numbers into $(9)$, we see that
$$a>Z^{-1}\times 111.4\,\mathrm{\mu m}.
\tag{10}$$
To make calculations more practical, we want possibly smaller value of $a$, so that our eigenfunctions would have fewer oscillations. Let's then take Einsteinium$^\ddagger$ as the nucleus, so that $Z=99$, and let's take the smallest value we can (a fraction of percent smaller, for nicer numbers):
$$a=1.12\,\mathrm{\mu m}.
\tag{11}$$
Let's define some more values now. We'll choose semi-minor axis $b$ of the orbit to have a reasonable value — so that it's obviously non-circular, but not rendering pericenter too small:
$$b=\frac34a=0.84\,\mathrm{\mu m}.
\tag{12}$$
Then we have average angular momentum $M$ defined by
$$M=b\sqrt{-2mE}=1.145\times10^{-31}\,\mathrm{J\cdot s}=1086\hbar,
\tag{13}$$
where $E$ is total energy (which is negative, since we're interested in a bound state). It is given by
$$E=\frac{-\alpha}{2a}=-1.02\times10^{-20}\,\mathrm{J}=-0.063\,\mathrm{eV}.
\tag{14}$$
From the equality
$$E_n=-\frac{Z^2me^4}{8(2\pi\hbar)^2\varepsilon_0^2n^2}
\tag{15}$$
we can see that average value of principal quantum number $n$ is
$$\langle n\rangle=1447.5.
\tag{16}$$
We'll set initial average distance of the electron from the nucleus to the orbit's apocenter:
$$r_0=a\left(1+\sqrt{1+\frac{2EM^2}{m\alpha^2}}\right)=1.86\,\mathrm{\mu m}.
\tag{17}$$
For ease of calculation of projections of our initial wave packet on eigenstates, we'll represent it as a product of three functions:
$$\psi_0(r,\theta,\phi)=R(r)\Theta(\theta)\Phi(\phi),
\tag{18}$$
where these functions are defined as
$$\begin{align}
\Phi(\phi)&=\exp\left(-\frac{(r_0(\phi-\pi))^2}{4\sigma_0^2}\right)\exp\left(i\frac M\hbar \phi\right),
\tag{19}\\
\Theta(\theta)&=\exp\left(-\frac{(r_0(\theta-\pi/2))^2}{4\sigma_0^2}\right),
\tag{20}\\
R(r)&=\exp\left(-\frac{(r-r_0)^2}{4\sigma_0^2}\right).
\tag{21}
\end{align}$$
With these definitions, the orbit and the $3\sigma_0$ circle of the initial wave packet will look like:
Practical calculations
In principle, the above listed parameters are enough to construct the wavepacket we want. But if we don't have a super-fast computer, we may want to do some simplifications to the calculations.
For $\Phi$ it is trivial to find projections $\Phi_m$ on eigenstates: it's just a Fourier transform, and since our wave packet is very localized in angular space, we may even integrate over $\phi\in(-\infty,\infty)$ for simplicity, getting results very close to integration over $\phi\in[0,2\pi]$. The result will thus be
$$\Phi_m=
\int\limits_0^{2\pi} \Phi(\phi)e^{-im\phi}\,\mathrm d\phi\approx
\int\limits_{-\infty}^\infty \Phi(\phi)e^{-im\phi}\,\mathrm d\phi=\\
=\frac{2\sqrt\pi\sigma_0}{r_0}\exp\left(\frac{(m-M/\hbar)((M/\hbar-m)\sigma_0^2-i\pi r_0^2)}{r_0^2}\right).
\tag{22}$$
For our chosen values it'll suffice to use the projections with $m=1000,...,1170$, so that minimum value of the projection is $\sim10^{-4}$ (maximum is $\sim0.1$).
For $\Theta$ finding projections is not that trivial, and I had to resort to numerical integration. The integral is
$$\Theta_{l,m}=\int\limits_0^\pi \Theta(\theta)Y_l^m(\theta,0)\sin(\theta)\,\mathrm d\theta,
\tag{23}$$
where $Y_l^m(\theta,\phi)$ is the spherical harmonic function. Luckily, due to the nature of high-angular momentum spherical harmonics, we only have to use $l=m,...,(m+7)$. Additionally, when $l+m$ is odd, the integral vanishes, so this eases the calculation a bit.
For $R$ finding projections directly is the hardest problem numerically. The integral is
$$R_{n,l}=\int\limits_0^\infty R(r)F_{n,l}(r)r^2\,\mathrm dr,
\tag{24}$$
where $F_{n,l}$ is radial eigenfunction of a hydrogen-like system:
$$F_{n,l}(\rho)=\left(\frac{a_B}Z\right)^{-3/2}\sqrt{\frac{(n-l-1)!}{(n+l)!}}e^{-\frac rn}\left(\frac{2r}n\right)^l\frac2{n^2}L_{n-l-1}^{2l+1}\left(\frac{2r}n\right),
\tag{25}$$
where we denoted $r=\frac{Z\rho}{a_B}$, and $a_B$ is the Bohr radius. Laguerre function $L_n^a$ is given in the notation used in GNU Scientific Library and in Wolfram Mathematica.
Despite the hardness of calculations of $R_{n,l}$, we can satisfactorily fit it using the following formula for the case of $l=m_{\mathrm{min}}=1000$:
$$R_{n,1000}\approx(-1)^{n+l+1}\exp\left(-154.84 - 0.92643 n + 0.00149354 n^2 - 5.46295\times 10^{-7} n^3\right).
\tag{26}$$
Maximum of $R_{n,l}$ as a function of $n$ roughly coincides with the case when $F_{n,l}(r)$ has a farther turning point at $r=r_0$. Thus, to find $R_{n,l}$ for other values of $l$, we can simply move its maximum by the difference between values of $n$ corresponding to the farther turning points being equal to $r_0$ for the eigenstates $F_{n,l}(r)$. The turning points can be found from the equation
$$-\frac\alpha {r_0}+\frac{\hbar^2l(l+1)}{2mr_0^2}=-\frac{Z^2me^4}{8(2\pi\hbar)^2\varepsilon_0^2n^2}.
\tag{27}$$
Calculations show that for satisfactory precision we can take $n=1300,...,1600$.
Results
After the calculations described above we get a finite expansion of our wave packet in hydrogen-like eigenfunctions:
$$\operatorname{packet}(r,\theta,\phi,t)=\sum\limits_{m=1000}^{1170}e^{im\phi}\Phi_m\sum\limits_{l=m}^{m+7}Y_{l,m}(\theta,0)\Theta_{l,m}\sum\limits_{n=1300}^{1600}R_{n,l}F_{n,l}(r)\exp\left(-i\frac E\hbar t\right).
\tag{28}$$
Here's how the wave packet evolution looks for the slice of $\theta=\frac\pi2$ for $t\in[0,0.95T]$ (a longer version can be downloaded from here, 50.3 MiB):
We can see that our estimations were a bit off: the wave packet doesn't appear to spread in the radial direction as much. But it's in fact good news.
Another thing to note is very pronounced spreading of the wave packet due to uncertainty in radial position at apocenter, which becomes higher the further we go in time (visible in the longer animation). But this is a completely classical effect. If we put several classical particles along the $x$ axis, with the same angular momentum, mass and charge as for our electron, and let them move (neglecting interaction between them), we'll see the same spreading of their positions. Here's an animation of ten such particles superimposed on the wave packet:
At the end of the 4th orbital period (see the longest animation linked above) the packet spreads so much that its "head" begins to interfere with its "tail", and further in time it begins to constantly self-interfere.
$^\dagger$For derivations of formulas describing orbital motion see e.g. Landau & Lifshitz, "Mechanics", $\S15$ Kepler's problem.
$^\ddagger$Because it's the heaviest nucleus obtainable in macroscopic-scale quantities.
