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. You can use either (7) or (8) to determine the corresponding initial value of $\theta$, i.e., where along the ellipse you’re starting. The major axis of the ellipse is in the direction where $\theta=0$.