Арифметические действия над матрицами - Pascal ABC

Формулировка задачи:

Помогите, пожалуйста, исправить программу. Она должна производить действия над матрицами по такой формуле: С(конечная матрица)=А*(А²-B)-2(B+A)B Матрица А= (2 3 1) (-1 2 4) (5 3 0) Матрица В= (2 7 13) (-1 0 5) (5 13 21)

program lab7;
uses crt;
var i,j,k: integer;
a:array[1..3,1..3] of integer;
b: array [1..3,1..3] of integer;
c: array [1..3,1..3] of integer;
cs: array [1..3,1..3] of integer;
cr: array [1..3,1..3] of integer;
begin
b[1,1]:=2;b[2,1]:=-1;b[3,1]:=5;b[1,2]:=7;b[2,2]:=0;b[3,2]:=13;b[1,3]:=13; b[2,3]:=5;b[3,3]:=21;
a[1,1]:=2;a[2,1]:=-1;a[3,1]:=5;a[1,2]:=3;a[2,2]:=2;a[3,2]:=3;a[1,3]:=1; a[2,3]:=4;a[3,3]:=0;
write('A:');
for i:=1 to 3 do
 begin
writeln();
for j:=1 to 3 do
  begin
write(a[i,j],' ');
 end;
  end;
writeln();
write('B:');
for i:=1 to 3 do
 begin
writeln();
for j:=1 to 3 do
  begin
write(b[i,j],' ');
 end;
 end;
For i:=1 to 3 do
begin
for j:=1 to 3 do
cs[i,j]:=a[i,j]*(a[i,j]*a[i,j]-b[i,j]);
end;
For i:=1 to 3 do
begin
for j:=1 to 3 do
cr[i,j]:=2*(b[i,j]+a[i,j])*b[i,j];
end;
For i:=1 to 3 do
begin
for j:=1 to 3 do
 begin
 c[i,j]:=0;
 for k:=1 to 3 do
  c[i,j] := c[i,j] + cs[i,k] - cr[k,j];
 end;
 end;
writeln();
write('C:');
for i:=1 to 3 do
 begin
writeln();
for j:=1 to 3 do
  begin
write(c[i,j],' ');
 end;
  end;
writeln();
end.

В ответе должна получиться такая матрица С= (-75 -357 -788) (-40 -156 -359) (-147 -607 -1291) А у меня она не получается.. Что не так? Помогите исправить.

Решение задачи: «Арифметические действия над матрицами»

Листинг программы

program lab7;
uses crt;
const n=3;
type mtr=array[1..n,1..n] of integer;
const a:mtr=((2,3,1),(-1,2,4),(5,3,0));
      b:mtr=((2,7,13),(-1,0,5),(5,13,21));
var i,j,k: integer;
    a2,a3,a4,b2,b3,c:mtr;
begin
writeln('Матрица A:');
for i:=1 to n do
 begin
  for j:=1 to n do
  write(a[i,j]:4);
  writeln;
 end;
writeln('Матрица B:');
for i:=1 to n do
 begin
  for j:=1 to n do
  write(b[i,j]:4);
  writeln;
 end;
//A^2
for i:=1 to n do
for j:=1 to n do
 begin
  a2[i,j]:=0;
  for k:=1 to n do
  a2[i,j]:=a2[i,j]+a[i,k]*a[k,j];
 end;
//A^2-B
for i:=1 to n do
for j:=1 to n do
a3[i,j]:=a2[i,j]-b[i,j];
//A*(A^2-B)
for i:=1 to n do
for j:=1 to n do
 begin
  a4[i,j]:=0;
  for k:=1 to n do
  a4[i,j]:=a4[i,j]+a[i,k]*a3[k,j];
 end;
//2*(B+A)
for i:=1 to n do
for j:=1 to n do
b2[i,j]:=2*(b[i,j]+a[i,j]);
//2*(B+A)*B
for i:=1 to n do
for j:=1 to n do
 begin
  b3[i,j]:=0;
  for k:=1 to n do
  b3[i,j]:=b3[i,j]+b2[i,k]*b[k,j];
 end;
//C=.....
for i:=1 to n do
for j:=1 to n do
c[i,j]:=a4[i,j]-b3[i,j];
writeln('Матрица C=A*(A^2-B)-2(B+A)B');
for i:=1 to n do
 begin
  for j:=1 to n do
  write(c[i,j]:6);
  writeln
 end
end.

Объяснение кода листинга программы

Объявлены константы: n=3, a, b, i, j, k, a2, a3, a4, b2, b3, c.
Выводится матрица A.
Выводится матрица B.
Вычисляется матрица a^2.
Вычисляется матрица a^2-b.
Вычисляется матрица a*(a^2-b).
Вычисляется матрица 2*(b+a).
Вычисляется матрица 2(b+a)b.
Вычисляется матрица c=a*(a^2-b)-2(b+a)b.
Выводится матрица c.