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))
$$
Short answer: Use the radial acceleration component $a_\mathrm r$ of the total $\bf a$ by the central force and then deduce the SHM equation using the properties of circular orbital motion (effective potential energy curves can also be used) that $\dot r= 0$ and keeping in mind the fact that SHM would only occur when the circular orbit is stable.
Here the approach is elaborately discussed:
The radial component of the central force $\bf F$ is given by
$$\mathbf F\cdot \mathbf e_\mathrm r~=~ m\left[\ddot r - r(\dot \theta)^2\right]\;.\tag 1 $$
Due to the small perturbation, the radius of the circular orbit changed from $r_0$ to $r(t) ~=~r_0 + \rho(t)~~_; ~~~\rho(0)\ll r_0\;.$
Using $\mathit l= mr^2\dot \theta=\textrm{const}$, $(1)$ becomes
$$\ddot\rho(t)- \frac{l^2}{m^2(r_0+ \rho(t))^3}~=~ \frac{F_\mathrm r(r_0+ \rho(t))}{m}~\tag {1-a} $$
Now, \begin{align}(r_0+ \rho(t))^{-3}&=r_0^{-3}\left(1+ \frac{\rho}{r_0}\right)^{-3}\\ &\approx r_0^{-3}\left(1-3\frac{\rho}{r_0}\right)\\ F_\mathrm r(r_0+\rho)&= F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho+ \ldots \\ &\approx F_\mathrm r(r_0) +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\rho\end{align}
So,
$$\ddot\rho- \frac{l^2}{m^2r_0^3}\left(1-3\frac{\rho}{r_0}\right)= \frac{F_\mathrm r(r_0)}m +\frac{\mathrm dF_\mathrm r}{\mathrm dr}\bigg|_{r=r_0}~\cdot \frac{\rho}m\tag{1-b} $$
From the orbit-equation$^\ddagger$, we get
$$\mathit l^2=- mr_0^3 F_\mathrm r(r_0)_; $$
So,
$$\ddot \rho+ \left[\frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\right]\rho~=~0\;.$$
This can be written as
$$\ddot \rho + \omega^2 \rho ~=~0\;.\tag 2$$
Now, this is a second-order differential equation for simple harmonic oscillator with frequency $\omega\;,$
where $$\omega^2 ~=~ \frac{3\mathit l^2}{m^2r_0^4}- \frac{1}{m}\frac{\mathrm dF_\mathrm r}{\mathrm dt}\bigg |_{r=r_o}\;.$$
Unless, the circular orbit is unstable, a simple harmonic radial oscillation about $r=r_0$ would ensue (i.e. SHM would occur if $\omega^2\gt 0\;.$)
$\ddagger$ We would derive the orbit-equation from $(1)$ expressing $r$ and its derivatives in terms of $\theta\;.$
\begin{align}\dot r &= \frac{\mathrm dr}{\mathrm d\theta}\frac{\mathrm d\theta}{\mathrm dt}\\ \ddot r&=\frac{\mathrm d}{\mathrm dt}\left(\frac{\mathrm dr}{\mathrm d\theta}\dot \theta\right)\\&=\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta \end{align}
Let
$$\mathbf H \stackrel{\text{def}}{=}
\mathbf r\times \dot{\mathbf r} =\frac{\mathbf L}m \;.$$
Therefore
\begin{align}\dot \theta &= \frac{|\mathbf H|}{r^2}\\ \implies \ddot \theta &=-\frac{2|\mathbf H|}{r^3} \, \dot r\;.\end{align}
Putting all these in $(1),$ we get
\begin{align}\frac{\mathrm d^2r}{\mathrm d\theta^2}~\dot\theta ^2 + \frac{\mathrm dr}{\mathrm d\theta}~\ddot \theta- r(\dot \theta)^2& = \frac{F_\mathrm r}{m} \\ \implies\frac{\mathrm d^2r}{\mathrm d\theta^2}~\left(\frac{|\mathbf H|}{r^2}\right)^2 +\frac{\mathrm dr}{\mathrm d\theta}~\left(-\frac{2|\mathbf H|}{r^3}\dot r\right) - r ~\left(\frac{|\mathbf H|}{r^2}\right)^2 &= \frac{F_\mathrm r}{m}\\ \implies \left(\frac{|\mathbf H|}{r^2}\right)^2\left[\frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r\right]&=\frac{\mathrm F(r)}{m}\\ \implies \frac{\mathrm d^2r}{\mathrm d\theta^2} - 2\frac{\left(\frac{\mathrm dr}{\mathrm d\theta}\right)^2}{r}-r&= \frac{r^4 F_\mathrm r}{|\mathbf H|^2 m}\;.\tag{i}\end{align}
$(\rm i)$ is the orbital equation.
Now, applying the facts that for circular orbits, $\dot r,~\ddot r~=~0\,,$ we get the value of $$|\mathbf H|~=~\sqrt{ -\frac{r^3F_\mathrm r}{m}}\;,$$ where $r~=~r_0\;.$
Best Answer
For the case of hyperbolic orbit, you use the hyperbolic trigonometric function. So your equation for Hyperbolic Mean anomaly becomes:
$$M_h = e\sinh F - F$$
where $F$ is your hyperbolic eccentric anomaly, which is analogous to eccentric anomaly for ellipse and $M_h$ is still given by the same expression. The angle measurement for hyperbolic trigonometric identities are not the same but analogous, so using the center of hyperbola which exists between the real hyperbola and it's image, you get the sinh function value. Assuming the object to be at distance $r$ from the primary focus and center of hyperbola being the reference origin $\sinh F$ will be given by:
$$\sinh F = \dfrac{y}{b}$$
where $y$ is the perpendicular distance from the focal axis of the object on trajectory, and $b$ is the semi-minor axis of the conic.
As for the semi-major axis it can be obtained by orbit equation at periapse for hyperbola which is $a = \dfrac{h^2}{\mu} \dfrac{1}{e^2-1}$, which is positive.
For parabolic orbit you do not define Mean Anomaly, because of rightly mentioned absence of geometric center, instead use Baker's equation $2\sqrt{\dfrac{\mu}{p^3}}\left(t-\tau\right) = \tan\dfrac{\theta}{2} +\dfrac{1}{3}\tan^3\dfrac{\theta}{2}$ where $\tau$ is your periapse passage time and $t$ is the elapsed time since periapse, and $p$ is the semi-latus rectum which is well defined for a parabola. There are solutions available to this equation that gives the position of the object at any given time given some input parameters.
Hope this helps.