Yes you can do that (easily) with Wolfram Language (aka Mathematica.)

Here is an example:

    mat1 = Table[
       RandomChoice[{
           Quantity[RandomReal[], "Meters"],
           Quantity[RandomReal[], "Seconds"],
           Quantity[RandomReal[], "Meters/Seconds"],
           Quantity[RandomReal[], "Meters/Seconds^2"]
         }] , 3, 3];
    mat1 // MatrixForm
    
    mat1 . Transpose[mat1]

I whipped up this in raku:

  use Physics::Measure :ALL;

  sub infix:<·>(@x1, @x2) {
    die "Incompatible dimensions."
            unless @x1 == @x2[0] && @x1[0] == @x2;

    [for ^@x1 -> $m {
        [for ^@x1 -> $n {
            [+] @x1[$m;*] >>*<< @x2[*;$n]
        }]
    }]
  }

  my @m = [[1m,2m],[3m,4m]]; 

  say @m · [Z] @m;     #[[5m^2 11m^2] [11m^2 25m^2]]