Integrating two ODE's in MATLAB

Question

I'm trying to integrate these two equations:

$$\frac{d}{da}(\rho a^3) = -3w\rho a^2$$

$$((\frac{da}{dt})/a)^2 + k/a^2 = \frac{8\pi G\rho}{3}$$

where $a==a(t)$, $\rho==\rho(a)$, $w=w(\rho)$.

For now, I am taking $w$ to be a constant - later I would like to express $w$ as a function of $\rho$. I found out $$\frac{d\rho}{dt} = \frac{d\rho}{da} \frac{da}{dt}$$ to get 'ode1' and expressed $\rho$ as a function of time $t$ instead. 'ode2' is just the second equation. So, my variable of integration is $t$.

Code:

---

syms  a(t) rho(t) w(rho)

G=1;
w=1/3;
k=0;

ode1 = diff(rho) == -3*rho*(w+1)*sqrt(8*pi*G/3*rho);
ode2 = (diff(a)/a)^2 + k/a^2 == 8*pi*G/3*rho;

odes = [ode1;ode2];
Sol = dsolve(odes)
aSol(t) = Sol.a;
rhoSol(t) = Sol.rho

---

Output:

Sol = 
  struct with fields:
      a: [1×1 sym]
    rho: [1×1 sym]


aSol(t) =
  C5


rhoSol(t) =
  0

K>> 

This is of course not the solution. Can someone please tell me where I'm going wrong. Thanks.

I am open to numerically integrating these as well