var
  a: array[1 .. 1000, 1 .. 2000] of real;
begin
  a[1, 1] := 0;
end.

Array[1 .. 1000, 1 .. 2000] Of Real

j := 1;
while not eof(f) do begin
  for i := 1 to {количество чисел в строке} do begin
    read(f, a[j]); inc(j);
  end;
  readln(f);
end;

procedure EnterMatrixFromFile (var Matr: Matrix);

const
n=5;
m=6;
var
  i,j: integer;
  f:text;
  TextElem:string;

 begin
 assign(f,'matrixdat.txt');
 reset(f);
 i := 1;
  while not eof(f) do
        begin
         while not eoln(f)do
            begin
               j:=1;
               read(f, Matr[i,j]);
               inc(j);
            end;
           readln(f);
           inc(i);
        end;

{$N+}
uses crt;

const n=3;

type

   TArr  = array[1..n,1..n] of Extended;
   TFile = Text;

var
   f_: TFile;
   x_: TArr;

procedure ReadFromFile(var x: TArr; var f: TFile);
var

  e,i,j,t,back  : integer;
          tempE : extended;
        s,tempS : string;

begin

   i:=1; j:=0;

   while not(eof(f)) do begin
      readln(f, s);
      t:=1;
      while(t<=length(s)) do begin
          while(s[t]=#32) do
           inc(t);

          if t<=length(s) then begin
              back:=t;
              while(t<=length(s))and(s[t]<>#32) do
               inc(t);

              tempS:=copy(s, back, t-back);

              Val(tempS,tempE,e);

              if j=n then begin
                  j:=1;
                  inc(i);
              end
               else
                inc(j);

              x[i,j]:=tempE;

          end;

      end;

   end;

end;

ix := 1; iy := 1;
while not eof(f) do begin
  for i := 1 to {количество чисел в строке} do begin
    read(f, a[ix, iy]); inc(iy);
  end;
  readln(f); inc(ix)
end;

ix := 1; iy := 1;
while not eof(f) do begin
  while not eoln(f) do begin
    read(f, a[ix, iy]); inc(iy);
  end;
  readln(f); inc(ix)
end;

program ReadMatrixFromFile;
 const
  NMAX = 100;
  MMAX = 100;
type
  Matrix = array [1..NMAX, 1..MMAX] of real;
var
  M: Matrix;
  

(***************************************************
            Процедура ввода матрицы из файла
****************************************************)
procedure EnterMatrixFromFile (Matr:Matrix);

var
  n,m,i,j: integer;
  f:text;


 begin
 assign(f,'matrixdat.txt');
 reset(f);
 i := 1;

 while not eof(f) do
   begin
   j := 1 ;
  while not eoln(f) do
         begin
             read(f,Matr[i, j]);
             writeln('Matr[',i,',',j,']',Matr[i,j]);
             inc(j);
         end;
       n:=i;
       m:=j-1;
       writeln('n=',n,',','m=',m);
       readln(f); inc(i);
   end;

   close(f);
 end;

 begin
        EnterMatrixFromFile(M);
        readln;
 end.

procedure EnterMatrixFromFile(Matr: Matrix; var size1, size2: integer);

Кстати, изменил размерность до 1000Х2000 программа перстала работать.

procedure Funky(Var Matrix_: Matrix; nn, mm: integer);

procedure PrintMatrix (Matrix_: Matrix; nn, mm: integer);

procedure PrintMatrix (Var Matrix_: Matrix; nn, mm: integer);

procedure PrintMatrixToFile (Matrix_: Matrix; nn: integer);

procedure PrintMatrixToFile (Var Matrix_: Matrix; nn: integer);

procedure Funky(Matrix_: Matrix; nn, mm: integer);
var
  i,j,k,p: integer;
  Sum1,Sum2,Sum3: real;
  Funk:Matrix;
...

var
  Funk:Matrix;

procedure Funky(Var Matrix_: Matrix; nn, mm: integer);
var
  i,j,k,p: integer;
  Sum1,Sum2,Sum3: real;
...

Sum3:=Sum3+((1/(Matrix_[k,j]+Matrix_[p,j]))*
  (Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2)*
  (Matrix_[k,j]/Sum1-Matrix_[p,j]/Sum2));

If (Matrix_[k,j]+Matrix_[p,j])<>0 then ...

Function Simps(b:extended;m:integer):extended;
var
  sum, h, a: extended;
  N, j:Integer;
begin
  N:=1000; a := 0; { так как у тебя интервал начинается с 0 }
  h:=(b - a)/N;
  sum:=0.5*(power(b,-1+m/2)*exp(-b/2));
  for j:=1 to N-1 do
    sum:=sum+(j mod 2+1)*power(j*h,-1+m/2)*exp(-j*h/2);
  Simps:=2*h*sum/3
end;

 const
  NMAX = 20;
  MMAX = 20;

function Gamma(x: extended): extended;
{ ... }
begin { gamma function }
  If x=0 then gamma := 1;
  If x>0.0 then begin { ... } end
  else { x<0 } begin
    { ... }
    gam:=gamma(y);{ recursive call }
    for i:=0 to j-1 do
      gam:=gam/(x+i); { <--- Здесь !!! }
    Gamma:=gam
  end { x<0 }
end; { gamma function }

If x = 0 Then Begin Gamma := 1; Exit End;
{ Далее - по тексту... }

function Gamma(x: extended): extended;
const
pi	= 3.1415926;

var
i,j	: integer;
y,gam	: extended;

begin		{ gamma function }
If x=0 then gamma:=1;
  if x>0.0 then
    begin
      y:=x+2.0;
      gam:=sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
      Gamma:=gam/(x*(x+1))
    end;
  if  x<0 then
    begin {x<0}
      j:=0;
      y:=x;
      repeat
	j:=j+1;
	y:=y+1.0
      until y>0.0;
      gam:=gamma(y);		{ recursive call }
      for i:=0 to j-1 do
	gam:=gam/(x+i);
      Gamma:=gam
    end; { x<0 }
end;		{ gamma function }