We study three computer algebra systems, namely SageMath (with SageManifolds package), Maxima (with ctensor package) and Python language (with GraviPy module), which allow tensor manipulation for general relativity calculations along with general algebraic calculations. We present a benchmark of these systems using simple examples. After the general analysis, we focus on the SageMath and SageManifolds system to derive, analyze and visualize the solutions of the massless Klein-Gordon equation and geodesic motion with Hamilton-Jacobi formalism. We compare our numerical result of the Klein-Gordon equation with the asymptotic form of the analytical solution to see that they agree. * Man = Manifold (4 , 'Man ' , r '\ mathcal { M } ') 4 5 # Define the parameter " M " ( mass ) : 6 M = var ( 'M ') 7 8 # Define the coordinates { t =0 , r =1 , theta (= th ) =2 , phi =3} with ranges 9 # ( BL = Boyer -Lidquist ) 10 BL . = Man . chart (r ' t r :(0 ,+ oo ) th :(0 , pi ) :\ theta ph :(0 ,2* pi ) :\ phi ') 11 12 # Define the metric " g " on manifold " Man ": 13 g = Man . lorentzian_metric ( 'g ') 14 15 # Enter the metric components : 16 g [0 ,0] = (1 -(2* M ) / r ) 17 g [1 ,1] = -1/(1 -(2* M ) / r ) 18 g [2 ,2] = -r^2 19 g [3 ,3] = -( r * sin ( th ) )^2 20 21 # Einstein tensor 40 ET = g . ricci () -(1/2) * g * g . ricci_scalar () 41 # Display all components 42 ET . set_name ( 'E ') 43 show ( ET . display () ) 44 # Display a single component 45 show ( ET [1 ,2]) ¦ ¥ The lines between 16-19 defines the Schwarzschild metric. In order to define the Kerr metric, we need to change this part with § ¤ 16 a = var ( 'a ') 17 rhosq = r^2+( a^2) * cos ( th )^2 18 Delta = r^2 -2* M * r + a^2 19 20 g [0 ,0] = (1 -(2* M * r ) / rhosq ) 21 g [1 ,1] = -rhosq / Delta 22 56 # We can analyze the Klein -Gordon equation 57 # to see how it is separated into 58 # radial and angular parts . 59 60 # Common factors will be divided 61 # This part should be given by the user 62 divideKGby = exp ( -I * omega * t ) * exp ( I * k * ph ) * sin ( th ) * R * S 63 finalKG = expand ( KG / divideKGby ) 64 65 # Extract the operands in the expression : 66 fkgops = finalKG . operands () 67 68 # Find radial and angular parts : 69 # lambda_aux is the separation constant 70 var ( ' lambda_aux ') 71 KGradialpart = lambda_aux 72 KGangularpart = -lambda_aux 73 for term in fkgops : 74 if diff ( term , r ) ==0: 75 KGangularpart = KGangularpart + term 76 else : 77 KGradialpart = KGradialpart + term 78 79 KGradialpart = expand ( KGradialpart * R) . simplify_full () . collect ( R ) . collect ( diff (R , r ) ) . collect ( diff (R ,r , r ) ) 80 KGangularpart = expand ( KGangularpart * S ) . simplify_full () . collect ( S ) . collect ( diff (S , th ) ) . collect ( diff (S , th , th ) ) 81 110 111 # Substitute numerical values for parameters 112 radial2 = radial2 . subs ( M =0.1 , omega =0.2 , k =2.0 , lambda_aux =2.0) 113 13 114 # Solve the system and plot the solution 115 radsol = desolve_system_rk4 ([ radial1 , radial2 ] ,[R , R_aux ] , ics =[0.3 ,1 ,0.5] , ivar ...