### Can't substitute H=dot(a)/a in SageMath

Hello all,
I'm working with the Friedman equations, and I've gotten them down to the common form presented in terms of the scale factor. MWE incoming:

```
M = Manifold(4, 'M', structure='Lorentzian')
fr.<t,r,th,ph> = M.chart(r't r:[0,+oo) th:[0,pi]:\theta ph:[0,2*pi):\phi')
var('G, Lambda, k', domain='real')
a = M.scalar_field(function('a')(t), name='a')
rho = M.scalar_field(function('rho')(t), name='rho')
p = M.scalar_field(function('p')(t), name='p')
g = M.metric()
g[0,0] = -1
g[1,1] = a*a/(1 - k*r^2)
g[2,2] = a*a*r^2
g[3,3] = a*a*(r*sin(th))^2
nabla = g.connection()
Ricci = nabla.ricci()
Ricci_scalar = g.ricci_scalar()
E1 = Ricci - Ricci_scalar/2*g + Lambda*g
print("First Friedmann equation:\n")
E1[0,0].expr().expand() == 0
```

With the result being

```
-Lambda + 3*diff(a(t), t)^2/a(t)^2 + 3*k/a(t)^2 == 0
```

Which is great, but I really want this to be in terms of the Hubble paramter $H=\dot{a}/a$. I can't make that substitution happen:

```
var('h')
eq1.substitute(diff(a.expr(),t)/a.expr()==h)
```

just spits out the same thing,

```
-Lambda + 3*diff(a(t), t)^2/a(t)^2 + 3*k/a(t)^2 == 0
```

Any ideas?