Решение системы нелинейных уравнений методом Ньютона

function F1(x,y:real):real;
begin
F1:=2-x*x-y*y;
end;
function F2(x,y:real):real;
begin
F2:=1-x*x-y*y;
end;

program MyGaus;
uses
crt;
const
nmax = 5; {кол-во уравнений}
type
t_matrix = array[1..nmax,1..nmax] of real;
t_answer = array[1..nmax] of real;
 
var
i,j,k,c,count   :byte;
arr  :t_matrix;
ans,b,x0     :t_answer;
tmp,tmp1,tmp2,E :real;
function a11(x1:Real): Real; {Частная производная по Х первого уравнения системы}
begin
a11 := -1/(2 * cos(x1) * sin(x1));
end;
function a21(x1:Real): Real; {Частная производная по Х второго уравнения системы}
begin
a21 := 200 * x1;
end;
function p1(x1,y1:Real): Real; {Вектор правой части для первого уравнения системы}
begin
p1 := sin(x1)/cos(x1) - y1;
end;
function P2(x1,y1:Real): Real; {Вектор правой части для второго уравнения системы}
begin
p2 := 100 * x1 * x1 - y1 - 10;
end;
begin
clrscr;
E := 0.0000001;
k:= 2; {кол-во уравнений в системе }
writeln('Vvedite priblizjenie ');
for i:=1 to k do begin
write('x',i,' = ');
read(x0[i]);
end;
count :=0;
repeat 
count :=count+1;
tmp1 := x0[1];
tmp2 := x0[2];
x0[1] := ans[1] + tmp1;
x0[2] := ans[2] + tmp2;
 
b[1] := p1(x0[1],x0[2]);
b[2] := p2(x0[1],x0[2]);
arr[1,1] := a11(x0[1]);
arr[2,1] := a21(x0[2]);
arr[1,2] := -1;
arr[2,2] := -1;
 
for i:=1 to k-1 do begin
for j:=i+1 to k do begin
tmp:= arr[i,i]/arr[j,i];
for c:=i to k do begin
arr[j,c] := arr[j,c] * tmp - arr[i,c];
end;
b[j] := b[j] * tmp - b[i];
end;
end;
ans[j]:=b[j]/arr[j,j];
{Oбратный ход - поиск корней}
for i:=k-1 downto 1 do begin
tmp:=b[i];
for j:=i+1 to k do begin
tmp:= tmp - ans[j] * arr[i,j];
end;
ans[i]:= tmp/arr[i,i];
end;
until (abs(x0[1] - tmp1) < E);
writeln('count = ',count);
for i:=1 to k do writeln('X',i,' = ',ans[i]);
readkey;
end.