Here is an algorithm you can implement.

Let the small orbiting body have mass $m$. You know its initial position and velocity vectors, $\mathbf{r}_0$ and $\mathbf{v}_0$, in an arbitrary Cartesian coordinate system.

I’ll assume that you want to treat the “significantly heavier” body with mass $M$ as being stationary and consider it to be the origin of your coordinate system. So in the formulas below I have assumed that $M\gg m$.

From the initial position and velocity you know the (constant) angular momentum,

$$\mathbf{L}=m\mathbf{r}_0\times\mathbf{v}_0\tag1.$$

This vector is perpendicular to the plane of the orbit, so now you know the orbital plane.

From the [*vis-viva* equation](https://en.wikipedia.org/wiki/Orbital_mechanics#Velocity)

$$v^2=GM\left(\frac2r-\frac1a\right)\tag2$$

you can use $\mathbf{r}_0$ and $\mathbf{v}_0$ to find the semimajor axis $a$ of the ellipse in this plane.

From the semimajor axis you can find the (constant) energy using

$$E=-\frac{GMm}{2a}\tag3.$$

From the energy and the angular momentum you can find the [orbital eccentricity](https://en.wikipedia.org/wiki/Orbital_eccentricity)

$$e=\sqrt{1+\frac{2EL^2}{G^2M^2m^3}}\tag4.$$

At this point you know the plane of the ellipse, the size of the ellipse, and the eccentricity of the ellipse. The remaining unknown is the orientation of the ellipse in the plane.

To find this, use the orbital equation

$$r=a\frac{1-e^2}{1+e\cos\theta}\tag5$$

where $\theta$ is an angular coordinate around the axis defined by $\mathbf L$, and the expression for the angular momentum in polar coordinates in the orbital plane,

$$L=mr^2\dot\theta\tag6.$$

These equations give the velocity components in terms of the angle around the ellipse as 

$$v_r=\dot{r}=\frac{L}{ma}\frac{e\sin\theta}{1-e^2}\tag7$$

and 

$$v_\theta=r\dot\theta=\frac{L}{ma}\frac{1+e\cos\theta}{1-e^2}\tag8.$$

You know one velocity, $\mathbf{v}_0$, and can find its $r$ and $\theta$ components. The value of $\theta$ that satisfies both (7) or (8) -- call it $\theta_0$ -- tells you where along the ellipse you’re starting. The major axis of the ellipse is in the direction where $\theta=0$.

ADDENDUM:

Here is a complete numerical example! Use units in which $GM=m=1$, and let the initial position be

$$\mathbf{r}_0=(1,2,3)$$

and the initial velocity be

$$\mathbf{v}_0=\left(\frac12,\frac13,\frac14\right).$$

One finds

$$\mathbf{L}=\left(-\frac12,\frac54,-\frac23\right),$$

$$E=\frac{427-144\sqrt{14}}{2016}\approx -0.0554557,$$

and

$$e=\sqrt{1-\frac{325(144\sqrt{14}-427)}{145152}}\approx 0.86564.$$

The initial radial velocity is

$$v_{r,0}=\mathbf{v}_0\cdot\frac{\mathbf{r}_0}{|\mathbf{r}_0|}=\frac{23}{12\sqrt{14}}\approx 0.512251$$

and the initial velocity along the $\hat\theta$ direction is

$$v_{\theta,0}=\sqrt{\mathbf{v}_0^2-v_{r,0}^2}=\frac{5}{12}\sqrt{\frac{13}{14}}\approx 0.40151.$$

Solving for $\theta_0$, one finds

$$\theta_0=\pi-\tan^{-1}\frac{115\sqrt{13\sqrt{395929+93600\sqrt{14}}}}{184679}\approx 117.277\text{ degrees}.$$

(I used *Mathematica*.)

To represent the orbit, we introduce some useful unit vectors.

A unit vector perpendicular to the orbit is

$$\hat{\mathbf{z}}=\frac{\mathbf{L}}{|\mathbf{L}|}=\frac{1}{\sqrt{13}}\left(-\frac65,3,-\frac85\right)\approx (-0.33282,0.83205,-0.44376).$$

A perpendicular unit vector pointing toward the initial position is

$$\hat{\mathbf{x}}=\frac{\mathbf{r}_0}{|\mathbf{r}_0|}=\frac{1}{\sqrt{14}}(1,2,3)\approx (0.267261,0.534522,0.801784).$$

A third perpendicular unit vector is

$$\hat{\mathbf{y}}=\hat{\mathbf{z}}\times\hat{\mathbf{x}}=\frac{1}{\sqrt{182}}\left(\frac{61}{5},2,-\frac{27}{5}\right)\approx (0.904324,0.14825,-0.400275).$$

We rotate around $\hat{\mathbf{z}}$ by $\theta_0$ to make new unit vectors where $\hat{\mathbf{x}}'$ points along the major axis:

$$\hat{\mathbf{x}}'=\hat{\mathbf{x}}\cos{\theta_0}-\hat{\mathbf{y}}\sin{\theta_0}\approx (-0.926249,-0.376731,-0.0116847),$$

$$\hat{\mathbf{y}}'=\hat{\mathbf{x}}\sin{\theta_0}+\hat{\mathbf{y}}\cos{\theta_0}\approx (-0.176901,0.407143,0.896019).$$

The orbit is then

$$\begin{align}
\mathbf{r}&=r(\hat{\mathbf{x}}'\cos\theta+\hat{\mathbf{y}}'\sin\theta)\\
&=a\frac{1-e^2}{1+e\cos\theta}(\hat{\mathbf{x}}'\cos\theta+\hat{\mathbf{y}}'\sin\theta)\\
&=\left(\frac{-2.09049\cos\theta-0.399255\sin\theta}{1+0.86584\cos\theta},\frac{-0.850262\cos\theta+0.9189\sin\theta}{1+0.86584\cos\theta},\frac{-0.0263717\cos\theta+2.02238\sin\theta}{1+0.86584\cos\theta}\right)
\end{align}.$$

I leave it to you to verify, as I did, that this satisfies the initial conditions when $\theta=\theta_0$.

As further sanity checks, one finds that when $\theta=0$ the position is $(-1.1204,-0.455699,-0.0141339)$ and when $\theta=\pi$ it is $(15.5821,6.33768,0.196569)$. The former is periapsis, at a distance of $1.20961$, and the latter is apoapsis, at a distance of $16.8228$. These sum to twice the semimajor axis, $2a\approx 18.0324$, and are individually $a(1-e)$ and $a(1+e)$ with $e\approx 0.86564$.

Note: All of this gets you the correct ellipse in 3D space, parameterized by the angle-around-the-orbit $\theta$. It doesn’t tell you where the object is at a given *time*, which adds more complications because $\theta$ isn’t proportional to $t$. For information about the time-dependence, see the Wikipedia article “[Kepler’s equation](https://en.wikipedia.org/wiki/Kepler%27s_equation)”.