J: matrix([1, a, b, c, e, d, f, gg, x], [a, c, e, b, f, gg, d, h, j], [b, e, d, f, gg, x, h
, j, k], [c, b, f, e, d, h, gg, x, u], [e, f, gg, d, h, j, x, u, v], [d, gg, x, h, j, k, u, v, w]
, [f, d, h, gg, x, u, j, k, r], [gg, h, j, x, u, v, k, r, s], [x, j, k, u, v, w, r, s, t] );

W:transpose([h,x,u,j,k,r,v,w])$

doallmxops: false$
doscmxops: false$
detout: true$

inv:invert(J)$
numInv:num(inv)$
denomInv:denom(inv)$


rho_p:submatrix(9,numInv,9)$
epsilon_p:numInv[9,9]$
V_p:submatrix(9,col(numInv,9))$

/* Verify the inverse */

for i:1 thru length(J) do
		for k:1 thru length(J) do
			print(ratsimp(row(J,k).col(numInv/denomInv,i)));


rhopW:rho_p.W$
rhopWW:rhopW.W$
VpW2:expand(expand((transpose(V_p).W))^2)$
epsilon_p:numInv[9,9]$
ep:epsilon_p$

rhopWW_Cs1:ratcoef(rhopWW, s, 1)$
rhopWW_Cs2:ratcoef(rhopWW, s, 2)$
rhopWW_Ct1:ratcoef(rhopWW, t, 1)$
rhopWW_const:expand(rhopWW-rhopWW_Cs1*s-rhopWW_Cs2*s^2-rhopWW_Ct1*t)$


D:denomInv$
D_Ct1:ratcoef(D,t,1)$
D_Cs1:ratcoef(D,s,1)$
D_Cs2:ratcoef(D,s,2)$
D_const:expand(D-D_Cs1*s-D_Cs2*s*s-D_Ct1*t)$

VpW:transpose(V_p).W$
VpW_Cs1:ratcoef(VpW,s,1)$

/* ratcoef(VpW,s,2), ratcoef(VpW,t,1), ratcoef(VpW,t,2) are all 0 */

VpW_const:expand(VpW-VpW_Cs1*s)$

save("VpWAndCoefficients.maxima",VpW, VpW_Cs1,VpW_const);




/********************* Calculate f0, f1, f2, f3, and save to disk **************/
                                                                         

(%i1) loadfile("VpWAndCoefficients.maxima");                                                    
(%o1)                      VpWAndCoefficients.maxima                                            
(%i2) loadfile("rhopWW.maxima");                                                                
(%o2)                            rhopWW.maxima        
(%i3) loadfile("VpWAndCoefficients.maxima");
(%o3)                      VpWAndCoefficients.maxima
(%i4) loadfile("ep.maxima");
(%o4)                              ep.maxima                                          
(%i5) values;                                                                                   
(%o5) [D, D_const, D_Cs1, D_Cs2, D_Ct1, ep, rhopWW, rhopWW_const, rhopWW_Cs1,
                               rhopWW_Cs2, rhopWW_Ct1, VpW, VpW_Cs1, VpW_const]

(%i6) f0:(rhopWW_const-(1/ep)*(VpW_const)^2)/D_const$                                           
(%i7) f1:(rhopWW_Cs1-(2/ep)*(VpW_const*VpW_Cs1))/D_Cs1$                                         
(%i8) f2:(rhopWW_Cs2-(1/ep)*(VpW_Cs1*VpW_Cs1))/D_Cs2$                                           
(%i9) f3:rhopWW_Ct1/D_Ct1$                                                                      
(%i10) save("f0.maxima", f0);                                                          
(%i11) save("f1.maxima", f1);                                                          
(%i12) save("f2.maxima", f2);                                                          
(%i13) save("f3.maxima", f3);                                                          



/*********** Verification of Equation 2.9 ********************/
    
(%i1) maxima_tempdir;                                              
(%o1)                           /home/easwaran                     
(%i2) maxima_tempdir:"/tmp";                                       
(%o2)                                /tmp                          
(%i3) maxima_tempdir;                                              
(%o3)                                /tmp                          
(%i4) showtime:true$                                               
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.     
(%i5) loadfile("f0.maxima")$                                       
Evaluation took 1.1450 seconds (1.5050 elapsed) using 98.784 MB.   
(%i6) loadfile("f3.maxima");                                       
Evaluation took 0.3170 seconds (0.3390 elapsed) using 19.708 MB.   
(%o6)                              f3.maxima                       


(%i15) ratsimped_f0:ratsimp(f0)$
Evaluation took 4362.8270 seconds (30767.7050 elapsed) using 311416.516 MB.

(%i16) save("ratsimped_f0",ratsimped_f0)$
Evaluation took 0.7730 seconds (1.1550 elapsed) using 37.483 MB.

(%i17) ratsimped_f3:ratsimp(f3)$
Evaluation took 0.8820 seconds (0.9540 elapsed) using 88.876 MB.

(%i18) save("ratsimped_f3",ratsimped_f3)$
Evaluation took 0.9260 seconds (0.9260 elapsed) using 37.482 MB.

(%i19) loadfile("f1.maxima")$
Evaluation took 0.7030 seconds (1.3930 elapsed) using 45.003 MB.

(%i20) ratsimped_f1:ratsimp(f1)$
Evaluation took 1525.1340 seconds (2847.1480 elapsed) using 125212.420 MB.

(%i21) save("ratsimped_f1",ratsimped_f1)$
Evaluation took 1.1390 seconds (1.5640 elapsed) using 37.481 MB.

(%i22) loadfile("f2.maxima")$
Evaluation took 0.2170 seconds (0.8790 elapsed) using 13.996 MB.

(%i23) ratsimped_f2:ratsimp(f2)$
Evaluation took 237.6640 seconds (241.9020 elapsed) using 21942.856 MB.

(%i24) save("ratsimped_f2",ratsimped_f2)$
Evaluation took 0.7310 seconds (0.7310 elapsed) using 37.483 MB.

(%i25) ratsimped_f0/ratsimped_f3;
Evaluation took 0.0190 seconds (0.0330 elapsed) using 4.000 KB.

(%o25)                                 1

(%i26) ratsimped_f1/ratsimped_f3;
Evaluation took 0.0220 seconds (0.0220 elapsed) using 0 bytes.

(%o26)                                 1

(%i27) ratsimped_f2/ratsimped_f3;
Evaluation took 0.0200 seconds (0.0200 elapsed) using 4.000 KB.

(%o27)                                 1