diff --git a/lib/matrix.gi b/lib/matrix.gi
index 9134c1840d..15292d66e8 100644
--- a/lib/matrix.gi
+++ b/lib/matrix.gi
@@ -4587,134 +4587,6 @@ local M, n, p, vars, slackVars, i, id, bestMove,
    return [point,value];
 end);
 
-
-# can do better for matrices and large exponents, preliminary improvement,
-# will be further improved (FL)
-##  InstallMethod( \^,
-##      "for matrices, use char. poly. for large exponents",
-##      [ IsMatrix, IsPosInt ],
-##  function(mat, n)
-##    local pol, indet;
-##    # generic method for small n, break even point probably a bit lower,
-##    # needs rethinking and some experiments.
-##    if n < 2^Length(mat) then
-##      return POW_OBJ_INT(mat, n);
-##    fi;
-##    pol := CharacteristicPolynomial(mat);
-##    indet := IndeterminateOfUnivariateRationalFunction(pol);
-##    # cost of this needs to be investigated
-##    pol := PowerMod(indet, n, pol);
-##    # now we are sure that we need at most Length(mat) matrix multiplications
-##    return Value(pol, mat);
-##  end);
-##
-# next iteration, conjugate matrix such that it is often very sparse
-# (a companion matrix), could still be improved, maybe with kernel functions
-# for compact matrices (FL)
-BindGlobal("POW_MAT_INT", function(mat, n)
-  local d, addb, trafo, value, t, ti, mm, pol, ind;
-  d := NrRows(mat);
-  # finding a better break even point probably also depends on q
-  if n < 2^QuoInt(3*d,4) or not IsField(DefaultFieldOfMatrix(mat)) then
-    return POW_OBJ_INT(mat, n);
-  fi;
-  # helper function to build up a semi-echelon basis
-  addb := function(seb, v)
-    local rows, pivots, len, vv, c, pos, i;
-    rows := seb.vectors;
-    pivots := seb.pivots;
-    len := Length(rows);
-    vv := ShallowCopy(v);
-    for i in [1..len] do
-      c := vv[pivots[i]];
-      if not IsZero(c) then
-        AddRowVector(vv, rows[i], -c);
-      fi;
-    od;
-    pos := PositionNonZero(vv);
-    if pos <= Length(vv) then
-      if not IsOne(vv[pos]) then
-        vv := vv/vv[pos];
-      fi;
-      Add(rows, vv);
-      Add(pivots, pos);
-      seb.heads[pos] := len + 1;
-      return true;
-    else
-      return false;
-    fi;
-  end;
-  # this returns a base change matrix such that t*m*t^-1 is block triangular
-  # with companion matrices along the diagonal
-  # (could/should? be improved to return t^-1, t*m*t^-1 and the
-  # characteristic polynomial of m at the same time)
-  trafo := function(m)
-    local id, b, t, r, a;
-    id := m^0;
-    b := rec(vectors := [], pivots := [], heads := []);
-    t := [];
-    # maybe better start with a random vector?
-    for a in id do
-      r := addb(b,a);
-      if r = true then
-        repeat
-          Add(t, a);
-          a := a*m;
-          r := addb(b,a);
-        until r <> true;
-      fi;
-    od;
-    t := Matrix(t, m);
-    return t;
-  end;
-  # compared to standard method, we avoid some zero or identity matrices
-  # and we multiply with mat from left to take advantage of sparseness of mat
-  value := function(pol, mat)
-    local f, c, i, val, j;
-    f := CoefficientsOfLaurentPolynomial(pol);
-    c := f[1];
-    i := Length(c);
-    if i = 0 then
-      return 0*mat;
-    fi;
-    if i = 1 then
-      val := POW_OBJ_INT(mat, f[2]);
-      return c[1] * val;
-    fi;
-    val := c[i] * mat;
-    if not IsMutable(val[1]) then
-      val := MutableCopyMatrix(val);
-    fi;
-    i := i-1;
-    for j in [1..NrRows(mat)] do
-      val[j,j] := val[j,j]+c[i];
-    od;
-    while 1 < i  do
-      val := mat * val;
-      i := i - 1;
-      for j in [1..NrRows(mat)] do
-        val[j,j] := val[j,j]+c[i];
-      od;
-    od;
-    if 0 <> f[2]  then
-      val := val * POW_OBJ_INT(mat, f[2]);
-    fi;
-    return val;
-  end;
-  t := trafo(mat);
-  ti := t^-1;
-  mm := t * mat * ti;
-  pol := CharacteristicPolynomial(mm);
-  ind := IndeterminateOfUnivariateRationalFunction(pol);
-  pol := PowerMod(ind, n, pol);
-  mm := value(pol, mm);
-  return ti * mm * t;
-end);
-
-InstallMethod( \^,
-    "for matrices, use char. poly. for large exponents",
-    [ IsMatrixOrMatrixObj, IsPosInt ], POW_MAT_INT );
-
 InstallGlobalFunction(RationalCanonicalFormTransform,function(mat)
 local cr,R,x,com,nf,matt,p,i,j,di,d,v;
   matt:=TransposedMat(mat);