with maxima

/* FLRW metric: (c) CC-BY-SA B. Roukema (2011) */
load(ctensor);
init_ctensor();
csetup();
4;
y;
[t,r,theta,phi];
1;
1;
   -c^2 ;   /* t */
   a^2 / (1 -k*r^2) ;  /* r */
   a^2 * r^2 ;    /* theta */
   a^2 * r^2 * (cos(theta))^2;  /* phi */

/* Enter functional dependencies with the DEPENDS function or 'N' if none */
/* e.g. depends([a,d],x); */
/* depends([a],t); */
depends([a],t);
y;

lg;  /* print the metric again (l= lower indices), use  ug  for the inverse */

einstein(true);  /* "mixed" means e.g. G_a^b  */

cdisplay(ein); /* print the Einstein tensor again */

ein[2,2] - ein[3,3];  /* both should equal p, but are they really the same? */

ein[2,2] - ein[4,4];  /* both should equal p, but are they really the same? */

/* define the stress-energy tensor */
T : matrix([-rho * c^2, 0,0,0], [0,p,0,0], [0,0,p,0], [0,0,0,p]);
depends([rho,p],t);

/* Friedmann equation */
ein[1,1] = 8* %pi * G/c^4 * T[1,1];  /* [1,1] component of G = 8 pi T */
-c^2/3 *lhs(%) = -c^2/3 * rhs(%), expand;    /* simplify */

/* acceleration equation */
/* [2,2] component of G = 8 pi T, and 
   subtract previous equation left- and right-hand sides, and 
   factorise */
c^2 * ein[2,2]  + lhs(%) =  c^2 * 8* %pi * G /c^4 * T[2,2] + rhs(%), factor;

/* simplify */
-lhs(%)/2 = -rhs(%)/2 ;