Из-за чего возникает ошибка? [pascal]

[supersuper2]

Пишу на Паскале 7.0(Pascal)
Ошибка деления на ноль.
Devision by zero.

Задание: Решить СЛАУ(систему линейных уравнений) методом Гауса по схеме единственного деления.
Система уравнений 4x4. В программе записал расширенную матрицу.

Вот само задание и описание метода подробное:
http://disk.tom.ru/cpdtvue

program q;
uses crt;
const i=4;

a:array[1..i,1..i+1] of real =
((2, 1, -5.1, 1, 2.7),
(5.4, 5.5, 4, -8.5, 21.9),
(5.3, -1, 1, 5.2, -3.9),
(1, 5.2, 2.5, -1, 9.9));
{type b = array[1..i, 1..i+1] of double;
type x = array[1..i, 1..i+1] of double;}
var j: integer;
b:array[1..i, 1..i+1] of real;
x:array[1..i] of real;
a1:array[1..i, 1..i+1] of real;
b1:array[1..i, 1..i+1] of real;
a2:array[1..i, 1..i+1] of real;
b2:array[1..i, 1..i+1] of real;
a3:array[1..i, 1..i+1] of real;
begin
for j:=1 to i+1 do begin
b[1, i+1]:= a[1,i+1]/a[1, 1];
a1[i,i+1]:= a[i,j] - a[i,1]*b[1,i+1];

b1[2,j]:= a1[2,j]/a1[2,2];
a2[i,i+1]:= a1[i,j] - a1[i,2]*b1[2,i+1];

b2[3, i+1]:= a2[3, i+1]/a2[3, 3];
a3[4, i+1]:= a2[4, i+1] - a2[4,3]*b2[3, i+1];  end;

for j:=1 to i do begin
x[4]:= a3[4,5]/a3[4,4];
x[3]:= b2[3,5] - b2[3,4]*x[4];
x[2]:= b1[2,5] - b1[2,3]*x[3] - b1[2,4]*x[4];
x[1]:= b[1,5] - b[1,2]*x[2] - b[1,3]*x[3] - b[1,4]*x[4];
writeln('x=', x[i]);
end;



readln;
end.

[Stilet]

a1[i,i+1]:= a[i,j] - a[i,1]*b[1,i+1];

Может тут надо a1[i,j+1]:= a[i,j] - a[i,1]*b[1,i+1];
Ты с индексами запутался(лась)