The Bel-Robinson tensor B_{ijhk} is a covariant rank 4 tensor defined in terms of the Weyl tensor W_{ijhk} on a 4-dimensional manifold by (see, for example, Penrose and Rindler Vol. 1)

The Bel-Robinson tensor is totally symmetric: B_{ijhk} = B_{jihk} = B_{hjik} = B_{kjhi}.

The Bel-Robinson tensor is trace-free: g^{ij} B_{ijhk} = 0.

If g is an Einstein metric, that is, R_{ij} = Lambda g_{ij} (where R_{ij} is the Ricci tensor for the metric g and Lambda is a constant), then the covariant divergence of Bel-Robinson vanishes: g^{il} nabla_l(B_{ijhk}) = 0.  Here nabla_l denotes the covariant derivative with respect to the Christoffel connection for g.

The keyword argument indexlist = ind allows the user to specify the index structure for the Bel-Robinson tensor.  For example, with indexlist = ["con", "con", "con", "con"], the contravariant form B^{ijhk} is returned.

This command is part of the DifferentialGeometry:-Tensor package, and so can be used in the form BelRobinson(...) only after executing the commands with(DifferentialGeometry); with(Tensor); in that order.  It can always be used in the long form DifferentialGeometry:-Tensor:-BelRobinson.