Unable to use substitute_function in SageManifolds

edit

I'm trying to do a function substitution into one of the Einstein Field equations, using tools from SageManifolds. I am doing nearly the same thing as the tutorial on the Lemaitre-Tolman equations (https://nbviewer.jupyter.org/github/e...). However, using substitute_function does not seem to achive the result. The MWE is below. The setup is the following:

M = Manifold(4, 'R^4')
coord.<t,x,z1,z2> = M.chart(r't x z1 z2') #name of the chart is coord.
a=function('a')(t)
g = M.metric('g')
g[0,0] = -1
g[1,1] = a^2
g[2,2] = a^2*sin(x)^2
g[3,3] = a^2*cos(x)^2
nab = g.connection()
ric = nab.ricci()
SCAL = g.ricci_scalar()
phi = M.scalar_field(function('phi')(*coord), name='phi')
dphi=nab(phi)
dphiu=dphi.up(g,0)
T1=dphi['_a']*dphi['_b'] 
T2=g.inverse()['^ab']*T1['_ab']
T3=dphi*dphi+1/2*g*T2

Obviously, on my sheet I run all these and they display what's expected - let me know if you think more info is needed here, and I'll attach the sheet as well (sorry, turns out I don't have the rep for that). So now I define the Einstein equations and try to make a trivial substitution:

var('G')
Lam=var('Lam', latex_name=r'\Lambda')
E=ric - SCAL/2*g + Lam*g - (8*pi*G)*T3
E[0,0].expr().expand() == 0

The result is

-12*pi*G*diff(phi(t, x, z1, z2), t)^2 + 4*pi*G*diff(phi(t, x, z1, z2), x)^2/a(t)^2 - Lam + 3*diff(a(t), t)^2/a(t)^2 + 4*pi*G*diff(phi(t, x, z1, z2), z1)^2/(a(t)^2*sin(x)^2) + 4*pi*G*diff(phi(t, x, z1, z2), z2)^2/(a(t)^2*cos(x)^2) + 3/a(t)^2 == 0

Cool, so what I really want to do is make some assumptions about the scalar field, but I can't even figure out how to substitute a different, equivalent function. My first try:

F1=M.scalar_field(function('F1')(*coord),name='F1')
E[0,0].expr().substitute_function(phi,F1)

Results in

(4*pi*G*cos(x)^2*diff(phi(t, x, z1, z2), z1)^2 + (4*pi*G*diff(phi(t, x, z1, z2), z2)^2 - (12*pi*G*a(t)^2*diff(phi(t, x, z1, z2), t)^2 - 4*pi*G*diff(phi(t, x, z1, z2), x)^2 + Lam*a(t)^2 - 3*diff(a(t), t)^2 - 3)*cos(x)^2)*sin(x)^2)/(a(t)^2*cos(x)^2*sin(x)^2)

e.g. exactly the same thing (specifically, phi is not replaced with F1). Here are bunch of other tries that also return the same thing:

E_expr=E[0,0].expr()==0
E_expr.substitute_function(phi, F1)
F2=function('F2')(*coord)
E[0,0].expr().substitute_function(phi,F2)
E_expr.substitute_function(phi, F2)

Anyone have any insight into what is happening?