Нашел формулу.

Экспонентом матрицы P называется наименьшее натуральное число t, такое, что все элементы матрицы Pt строго положительны. Обозначим его через exp P.

о вычислении и формуле
http://www.spbstu.ru/public/m_v/N_002/Sushkova/par_07.html

Ozzя, сам что-ли кропал? А у нас это вроде в FAQ-е есть
Как вычислить заданный многочлен от матрицы A

Угу... Вот и я про это:
Многочлен = 1*A^0 + 1*A^1 + (1/2)*A^2 + (1/6)*A^3 + ...

Для матрицы


Точность невелика, результаты, сопоставимые с матлаб, достигают только при кол-ве итераций, равных 13.

Результаты, полученные программой:


Результаты, полученные Matlab:



Кстати, volvo, в процедуре возведения матрицы в степень - ошибка. Должно быть
For i := 1 To pow-1 Do

{ Порядок матрицы }
Uses Crt;
Const
  size = 3;
  eps = 10E-6;

Type
  TMatrix = Array [1 .. size, 1 .. size] Of real;

Function Factorial(I:LongInt):LongInt;
var
  j,fact:LongInt;
begin
  if i=1 then
    fact:=1
  else
    begin
     fact:=1;
     for j:=2 to i do fact:=fact*j;
   end;
Factorial:=fact;
end;

(* Умножение матриц *)
Procedure matrixMult(Var m: TMatrix; a, b: TMatrix);
Var
  i, j, k: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    Begin
      m[i, j] := 0;
        For k := 1 To size Do
          m[i, j] := m[i, j] + a[i, k] * b[k, j]
    End;
End;

(* Возведение матрицы в степень *)
Procedure matrixPower(Var m: TMatrix; a: TMatrix; pow: Integer);
Var
  i, j: Integer;
  T: TMatrix;
Begin
  If pow = 0 Then
    Begin
      For i := 1 To size Do
      For j := 1 To size Do
        m[i, j] := Byte(i = j);
      Exit
    End;
    move(a, T, SizeOf(T));
    For i := 1 To pow-1 Do
      matrixMult(T, T, a);
      move(T, m, SizeOf(T))
End;

(* Сложение матриц *)
Procedure matrixAdd(Var m: TMatrix; a, b: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    m[i, j] := a[i, j] + b[i, j]
End;

(* Умножение матрицы на число *)
Procedure matrixScale(Var m: TMatrix; a: TMatrix; f: real);
Var
  i, j, k: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    m[i, j] := a[i, j] / f
End;

(* Печать матрицы *)
Procedure matrixPrint(a: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
    Begin
      For j := 1 To size Do
        write(a[i, j]:16:4);
      WriteLn
    End
End;

Procedure SingleMatrix(Var a: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
  if i=j then
    a[i, j]:=1
  else
   a[i, j]:=0;
End;

Var
  a,EMatrix,a2: TMatrix;
  T2, Res: TMatrix;
  i, j,n: Integer;
  ExpMatr:TMatrix;
  fact2:real;
begin
  clrscr;
  randomize;
  For i := 1 To size Do
  For j := 1 To size Do
    begin
    a2[i, j] :=2;
    a[i,j]:=a2[i,j];
    end;
  WriteLn('Исходная матрица');
  matrixPrint(a2);

  { Single E-Matrix here }
  SingleMatrix(EMatrix);
{  MatrixAdd(ExpMatr,EMatrix,ExpMatr);
  MatrixAdd(ExpMatr,A,ExpMatr);}

  n:=2;
  repeat
    For i := 1 To size Do
      For j := 1 To size Do
         a2[i,j]:=a[i,j];


    writeln('Matrica v stepeni ',n);
    MatrixPower(a2, a2, n);
    matrixPrint(a2);
    writeln;
    writeln('Matrica v stepeni,podelennay na fact');
    fact2:=Factorial(n);
    MatrixScale(a2,a2,Fact2);
    matrixPrint(a2);
    writeln;
    writeln('Matrica v stepeni podelennay na fact + suuma');
    MatrixAdd(ExpMatr,ExpMatr,A2);
    matrixPrint(ExpMatr);
    writeln;
    n:=n+1;
{    readkey;}
  until n>13;
  MatrixAdd(ExpMatr,EMatrix,ExpMatr);
  MatrixAdd(ExpMatr,A,ExpMatr);
  {abs(delta) <= eps;}

  writeln;
  clrscr;
  matrixPrint(ExpMatr);
  WriteLn('Factorial from ',n,' = ',Factorial(n));
  readln;
end.