REAL ph = 1.2;
[0:3,0:3]REAL a, t, tp;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  a[i,j] := t[i,j] := tp[i,j] := 0 OD OD;
a[0,0] := a[3,3] := cosh(ph); a[0,3] := a[3,0] := sinh(ph);
a[1,1] := a[2,2] := 1;
t[0,0] := 1; t[3,3] := -1;

FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
FOR k FROM 0 TO 3 DO FOR l FROM 0 TO 3 DO
  tp[j,l] +:= a[i,j]*a[k,l]*t[i,k] OD OD OD OD;
FOR i FROM 0 TO 3 DO FOR j FROM 0 TO 3 DO
  IF tp[i,j] /= 0 THEN
    print(("tp[",i,",",j,"] = ",tp[i,j],newline)) FI OD OD
# =>
tp[         +0,         +0] = +1.00000000000000e  +0
tp[         +3,         +3] = -1.00000000000000e  +0
#