Here is an algorithm you can implement.
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.
Your initial data are two vectors, $\mathbf{r}_0$ and $\mathbf{v}_0$, and the mass $m$ of the orbiting body.
From these 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
$$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
$$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! Take $G=M=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)$. TheseThe 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”.
