405.sagews

t,r,theta,phi = coords = var('t,r,theta,phi')
alpha = function('alpha', r)
beta = function('beta', r)
GM = var('GM')

metric = matrix([[-e^(2*alpha),0,0,0],[0,e^(2*beta),0,0],[0,0,r^2,0],[0,0,0,r^2*sin(theta)^2]])
inverseMetric = metric^(-1)
manifoldDimension = rank(metric)
def christoffel(lambda_, mu, nu):
    return sum([1/2*inverseMetric[lambda_][sigma]*(diff(metric[nu][sigma], coords[mu]) + diff(metric[sigma][mu], coords[nu]) - diff(metric[mu][nu], coords[sigma])) for sigma in range(manifoldDimension)])
christoffels = [[[christoffel(lambda_, mu, nu) for nu in range(manifoldDimension)] for mu in range(manifoldDimension)] for lambda_ in range(manifoldDimension)]
def riemann(rho, sigma, mu, nu):
    return diff(christoffels[rho][nu][sigma], coords[mu]) - diff(christoffels[rho][mu][sigma], coords[nu]) + sum([christoffels[rho][mu][lambda_]*christoffels[lambda_][nu][sigma] for lambda_ in range(manifoldDimension)]) - sum([christoffels[rho][nu][lambda_]*christoffels[lambda_][mu][sigma] for lambda_ in range(manifoldDimension)])
riemann(1,3,1,3)