Гаусс Опции

[marwell]

доброго времени суток
скорее всего проблема от банальной невнимательности, или просто что то упускаю. Но найти эту ошибку пытаюсь уже второй день,и безрезультатно. В общем передаю процедуре решения СЛАУ методом гаусса массив значений. Эта процедура заведомо правильная, написана не мной, найдена в интернете) Но у меня в коде почему то возникает деление 0/0. Хотя когда я ввожу ту же матрицу в исходную процедуру гаусса, все решается корректно. u1..u5 ввожу соответственно -1, -2, 1, 2, 3. Из этих u1..u5 собирается матрица вида
x0+x1*u1+x2*u1^2+x3*u1^3+x4*u1^4=f(u1);
x0+x1*u2+x2*u2^2+x3*u2^3+x4*u2^4=f(u2);
...
x0+x1*u5+x2*u5^2+x3*u5^3+x4*u5^4=f(u5);
вот ее я и передаю для решения методом гаусса.
буду очень благодарен совету

type
    Tmass=array[1..5] of real;
    Tmatrix=array[1..5,1..5] of real;
var
  Form1: TForm1;
  C1,C2,C3,C4,i,j:integer;
  a,b:integer;
  x,y1,y2,y3,y4,y5:integer;
  u1,u2,u3,u4,u5:real;
  s,sum:double;
  x1:Tmass;

implementation

{$R *.dfm}

function step(x:double;y:integer):double;
var c: double;
begin
if y=0 then c:=1
else begin
c:=(-1)*exp(y*ln(abs(x)));
if round(y/2)=y/2 then c:=abs©; end;
step:=c;
end;

function f(x: double): double;
begin
  f :=C1*x*x*x+C2*x*x+C3*x+C4;
end;

procedure TForm1.Button1Click(Sender: TObject);
var x0,y0:integer;
  w, h: integer;
  x:double;
begin
...
C1:=StrToInt(Edit1.Text);
C2:=StrToInt(Edit2.Text);
C3:=StrToInt(Edit3.Text);
C4:=StrToInt(Edit9.Text);
u1:=StrToFloat(Edit6.Text);
u2:=StrToFloat(Edit7.Text);
u3:=StrToFloat(Edit8.Text);
u4:=StrToFloat(Edit10.Text);
u5:=StrToFloat(Edit11.Text);
...
end;

procedure TForm1.Button5Click(Sender: TObject);
procedure Gauss(ish_matr:TMatrix; free:TMass);
var
    n,k:integer;d:real;
begin
n:=5;
    for i:=1 to n-1 do
       begin
         for k:=i+1 to n do
           begin
             d:=ish_matr[k,i]/ish_matr[i,i]; //вот тут ошибка invalid floating point operation
             for j:=i+1 to n do
               ish_matr[k,j]:=ish_matr[k,j]-ish_matr[i,j]*d;
             free[k]:=free[k]-free[i]*d;
           end;
       end;
    for i:=n downto 1 do
       begin
         for j:=i+1 to n do
            free[i]:=free[i]-ish_matr[i,j]*x1[j];
         x1[i]:=free[i]/ish_matr[i,i];
       end;
       end;

var uz: array [1..5] of real;
    ish_matr: Tmatrix;{исходная матрица коэффициентов}
       free: Tmass;{свободные члены}
       n,k,v: integer;
       x2,z,r,g: double;
begin
n:=5;
uz[1]:=u1;
uz[2]:=u2;
uz[3]:=u3;
uz[4]:=u4;
uz[5]:=u5;
for i:=1 to n do begin
for j:=1 to n do begin
ish_matr[i,j]:=step(uz[i],j-1);
end;
end;
for i:=1 to n do
free[i]:=f(uz[i]);
Gauss(ish_matr,free);
....
end;
end.