As far as I can tell the true anomaly is the same type of angle used in the standard solution of the Kepler problem since there we assume the sun is at a foci.
When solving the equations of motion for a Keplerian orbit we obtain $r\left(\theta\right)=\frac{a\left(1-e^{2}\right)}{1\pm e\cos\theta}$ (- if $r\left(0\right)$ is through the origin and + if it is away from the origin) and we can express the velocity as $v=\sqrt{\mu\left(\frac{2}{r}-\frac{1}{a}\right)}$ where $\mu$ is the standard gravitational parameter which in the two body case is just $G\left(m_{1}+m_{2}\right)$.
The derivation of these formulae can be found in many mechanics textbooks. For example, Taylor - Classical Mechanics - chapter 8.
First of your equation of your semi-major axis seems correct to me.
Secondly there is no explicit equation of your position as a function of time. However there is one the other way around, so an equation for time as a function of your position:
$$
t(\theta_0,\theta_1)=\sqrt{\frac{a^3}{\mu}}\left[2\tan^{-1}\left(\frac{\sqrt{1-e^2}\tan{\frac{\theta}{2}}}{1+e}\right)-\frac{e\sqrt{1-e^2}\sin{\theta}}{1+e\cos{\theta}}\right]^{\theta_1}_{\theta_0}
$$
For small steps in time you could use a polynomial, since you can find any order time derivative of the true anomaly $\theta$:
$$
\omega(\theta)=\frac{\delta}{\delta t}\theta=\frac{\delta \theta / \delta \theta}{\delta t / \delta \theta}=\frac{1}{\delta t / \delta \theta}=\sqrt{\frac{\mu}{a^3\left(1-e^2\right)^3}}\left(1+e\cos{\theta}\right)^2,
$$
$$
\alpha(\theta)=\frac{\delta^2}{\delta t^2}\theta=\frac{\delta}{\delta t}\omega=\frac{\delta \omega / \delta \theta}{\delta t / \delta \theta}=\frac{-2e\mu\sin{\theta}}{a^3(1-e^2)^3}(1+e\cos{\theta})^3,
$$
$$
ect.
$$
So:
$$
\theta(t_0+\Delta t)=\theta_{t_0}+\Delta t\omega(\theta_{t_0})+\frac{1}{2}\Delta t^2\alpha(\theta_{t_0})+O(\Delta t^3)
$$
However for big steps in time this will not be an accurate method, but you could use it to construct an iterative numerical method by calculating the time from $\theta(t_0+\Delta t)$ to determine the new $\Delta t$:
$$
\theta_{n+1}=\theta_n+\Delta t_n\omega(\theta_n)+\frac{1}{2}\Delta t_n^2\alpha(\theta_n)
$$
$$
\Delta t_{n+1}=\Delta t_n-t(\theta_n, \theta_{n+1})
$$
In certain cases it might be that this method would be unstable (when $|\Delta t_{n+1}|>|\Delta t_n|$) in that case you could split up your time step into multiple steps, or increase the order of the polynomial. There are probably a lot of other numerical methods which are better/faster, but I am not an expert in this area, and this should work and give you in idea of how to solve this.
For closed orbits you might also be able to reduce the size of your time step by taking a modified modulo operation: $\Delta t_{new}=\Delta t-round\left[\frac{\Delta t}{T}\right]T$, where $T$ is equal to the period of the orbit.
Edit I
I was messing with this myself in MATLAB and noticed that I often encountered a problem. I think most of the time because the calculated time was complex because $\theta_{n+1}$ exceeded the maximum range, for example a hyperbolic orbit will never be able to reach $\theta=\pi$. This range is determined by the $\tan^{-1}$, so when $e\ge 1$ than $\theta_{max}=2\tan^{-1}\left(\frac{1+e}{\sqrt{e^2-1}}\right)$. Now it does seems to work reliable, unless I try to calculate something which is getting near the machine accuracy. By the way I do not know in which computer language you are programming this simulator, but there is a big change that it will not be able to handle complex numbers. So you would have to keep this in mind when performing calculations, since each calculation should return real numbers, but certain sub-results will imaginary.
Edit II
I also looked into other methods to solve this, such as Runge–Kutta, however I was not able to prevent unstable behavior (when the time step is to big). However a method which would always be stable is a sort of binary search algorithm, in which check if an angle $\theta_{n+1}=\theta_n + \Delta\theta$ would yield a step in time smaller than you desired step and if so save that value. And after each iteration the step $\Delta\theta$ will be divided by two. For this you do have to choose the initial $\Delta\theta$ correctly, namely $\Delta\theta_0=\tan^{-1}\left(\frac{1+e}{e^2-1}\right)$ when $e\ge 1$ and $\Delta\theta_0=\pi$ when $e<1$. The convergent rate is not as good as that of the previous methods, especially when $e$ is big and the desired time step will result in a very high radius. So you might want to switch methods there, but I do not know when the other method will be stable. Another option might be to use the first method, but to perform the calculations relative to the radius, $r$, instead $\theta$:
$$
t(r)=\sqrt{\frac{a^3}{\mu}}\left(2\tan^{-1}{\sqrt{\frac{r-a(1-e)}{a(1+e)-r}}}-\sqrt{e^2-\left(\frac{r}{a}-1\right)^2}\right)
$$
$$
\dot{r}=\sqrt{\frac{\mu}{a}}\frac{\sqrt{e^2a^2-(r-a)^2}}{r}
$$
$$
\ddot{r}=\mu\left(\frac{a(1-e^2)}{r^3}-\frac{1}{r^2}\right)
$$
$$
\dddot{r}=\sqrt{\frac{\mu^3}{a}}\frac{\sqrt{e^2a^2-(r-a)^2}}{r^5}(2r-3a(1-e^2))
$$
Best Answer
Obviously it is not as obvious as I thought!
The radial velocity curve is defined through 6 free parameters $$V_r(t) = K\left(\cos(\omega + \nu) +e \cos \omega \right) + \gamma,$$ where $K$ is the semi-amplitude, $\gamma$ is the centre of mass radial velocity, $\omega$ is the usual angle defining the argument of the pericentre measured from the ascending node and $\nu$ is the true anomlay, which is a function of time, the fiducial time of pericentre passage $\tau$, the orbital period $p$ and the eccentricity $e$.
To proceed you estimate what all these parameters are - i.e. an initial guess.
Then, for each time $t_i$ of a data point in your RV curve you:
Calculate the mean anomaly $$M(t) = \frac{2\pi}{p}(t - \tau),$$
Solve "Kepler's equation" $$M(t) = E(t) - e \sin E(t)$$ numerically to give $E(t_i)$, the eccentric anomaly.
Use $$\tan \frac{E(t)}{2} = \left(\frac{1+e}{1-e}\right)^{-1/2} \tan \frac{\nu(t)}{2}$$ to calculate the true anomaly $\nu(t_i)$.
Calculate $V_r(t_i)$
You then calculate some figure of merit (e.g. chi-squared) and go through an iterative process to adjust the parameters and optimise the fit of model to data.
A more sophisticated discussion can be found in this paper by Beauge et al.