Numerically solving for cosmological scale factor in mathematica I would like to solve the scale factor $a(t)$ in a FRW metric.
The Friedmann equation for the scale factor is:
$H^2 \equiv \left( \frac{\dot{a}}{a} \right)^2 = 
H_0^2\left[\Omega_R \left(\frac{a(t)}{a(t_0)}\right)^{-4}+\Omega_M \left(\frac{a(t)}{a(t_0)}\right)^{-3} +\Omega_\Lambda  \right]\\ $
$
\frac{d^2a(t)}{dt^2} =   H_0^2\left[ 
\Omega_R
\frac{a(t_0)^4}{a(t)^2}
+ {\Omega_M}
\frac{a(t_0)^3}{a(t)}
+\Omega_\Lambda \right]$
I believe $a(t_0\equiv t=0)=1$ which represents the Universe today. $H_0 \simeq 1.44\times 10^{-42} \ { \text{GeV}} $ 
is the expansion rate measured today,
$\Omega_R\simeq 9.2\times 10^{-5}$ for radiation, 
$\Omega_M \simeq 0.31$ for matter, 
and  $\Omega_\Lambda\simeq 0.69$ for dark energy. In the above equation I take $a(t_0) = 1$ for the reason I outline above.
I am trying to solve this numerically in mathematica but I am unsure what to do for the boundary conditions as I am getting an error when I solve the above double differential equation. I take $a(0)=1$ i.e. the Universe today but what do I take for the other boundary condition. I am solving for earlier times so is this a negative $t$? i.e. is the other boundary condition $a(-100000)=0$. 
I suggest this as at the very start of the Universe $a(t_{\text{early}})= 0$. 
 A: I think you are making things complicated by writing the Friedmann equation in terms of of$\frac{a(t)}{a(t_0)}$. 
The general form is to write it in this way
$$\frac{H^2} {H_0^2} = \Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}$$
So this becomes 
$$H^2= H_0^2[\Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}]$$
$$\frac{\dot{a}(t)^2} {a(t)^2} =H_0^2[\Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}]$$
$$\dot{a}(t)^2=H_0^2[\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2]$$
let's take the square root
$$\dot{a}(t)=H_0 \sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}$$
From now on we can write
$$\frac{da}{dt}=H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}$$
Or 
$$\frac{da}{H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}} = dt$$
Let us call $$F(a) = \frac{1}{H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}}$$
So you have $$F(a)da = dt$$ which can be solved numerically.
From here you can easily get $t(a)$ (time as a function of the scale factor). But since we want $a(t)$ you just have to plot the data points as $x=t$ and $y=a$.
For the boundry conditions as many other poeple pointed out $a(t)$ should start from $0$. While doing the computation $0$ will cause an error since in the $F(a)$ we have $1/a^2$ etc. However you can avoid this by simply starting from $a = 0.0000...1$. For the upper limit you can put any number (non-negative) you want.

$H_0 \simeq 1.44\times 10^{-42} \ { \text{GeV}} $ 
  is the expansion rate measured today,

You should use  $H_0$ in terms of its normal units, $[H] = kms^{-1}Mpc^{-1}$. Or simply $H_0 \approx 70 km/s/Mpc$.
A: G. Smith is right on several points. 
One is that when you multiplied by $a(t)^2$ you forgot the third term in the right-hand-side of your second equation. It should be $a(t)^2 \Omega_\Lambda$
But the worst is the left-hand-side. 
It should be $(da/dt)^2$, not $d^2a/dt^2$. 
It is a first order equation (albeit of the second degree), not a second order equation. So you need only one boundary condition, namely $a(t)=1$ now, that is $a(t_0)=1$
You solve it by taking the positive square root of the right-hand side, and this is the appropriate value of $da/dt$.
How you solve this numerically, starting from $a(t_0)=1$ as single initial condition for a first order equation is your problem.
In fact, I did it recently, using the most primitive explicit integration scheme, forwards in time, but only for two or three billion years, to check a value I had found in the literature, and it checked perfectly. I forgot on which wiki page that value was, I am sorry to say. (Of course I completely neglected $\Omega_R$)
As for the units, $H_0$ has the dimension of the inverse of a time (velocity divided by distance). Therefore $H_0$ multiplied by the Planck constant would indeed be an energy. I did not check it, but a reasonable value for $H_0$ might be $1.44\ 10^{-42}$ GeV/$\hbar$, but you need to divide by ℏ to get a meaningful result.
Why on earth would anyone use a GeV/$\hbar$ as a unit for "inverse time" is hard to imagine. But at least it is meaningful. Just GeV  as unit of "inverse time" is meaningless..
