Решение системы нелинейных уравнений

#include <vcl.h>
#include <iostream.h>
#include <math.h>
#pragma hdrstop
 
#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
const int n=3;
int i,j;
double x,y,z,esp;
 
double f1(double x, double y, double z);
double f1(double x, double y, double z){
    return 0.1-pow(x,2)+2*y*z-x; }
double f2(double x, double y, double z);
double f2(double x, double y, double z){
    return -0.2+pow(y,2)-3*x*z-y;   }
double f3(double x, double y, double z);
double f3(double x, double y, double z){
    return 0.3-pow(z,2)-2*x*y-z;}
 
void __fastcall TForm1::Button2Click(TObject *Sender)
{
x=StrToFloat(Edit1->Text);
y=StrToFloat(Edit2->Text);
z=StrToFloat(Edit3->Text);
esp=StrToFloat(Edit4->Text);
 
double a[n][n], a1[n][n], b[n], det, invdet, f[n], K[n]; //a[3][3]-ìàòðèöà ßêîáè;
//a1[3][3]-îáðàòíàÿ ìàòðèöà; b[3]-ïðîèçâåäåíèå à1[3] íà f[3]; f[3]-ìàññèâ ôóíêöèé; K[3]-ìàññèâ ïåðåìåííûõ
double xk, yk, zk, xn, yn, zn, max; //èçìåíåíèå çíà÷åíèÿ ïåðåìåííûõ ïðè k=0,1,2..
int k=0; //êîë-âî èòåðàöèé
 
  //Ìàòðèöà ßêîáè
    {
        a[0][0]=(-2*x-1);
        a[0][1]=(2*z);
        a[0][2]=(2*y);
        a[1][0]=(-3*z);
        a[1][1]=(2*y-1);
        a[1][2]=(-3*x);
        a[2][0]=(-2*y);
        a[2][1]=(-2*x);
        a[2][2]=(-2*z-1);
 
     //Îïðåäåëèòåëü
     det=a[0][0]*a[1][1]*a[2][2]+a[0][1]*a[1][2]*a[2][0]+a[0][2]*a[1][0]*a[2][1]-
        a[0][2]*a[1][1]*a[2][0]-a[0][0]*a[1][2]*a[2][1]-a[0][1]*a[1][0]*a[2][2];
     invdet=1/det;
     //Îáðàòíàÿ ìàòðèöà
        for (i=0; i<n; i++){
         for (j=0; j<n; j++){
          a1[0][0]=a[1][1]*a[2][2]-a[1][2]*a[2][1]*invdet;
          a1[0][1]=-a[1][0]*a[2][2]+a[2][0]*a[1][2]*invdet;
          a1[0][2]=a[1][0]*a[2][1]-a[1][1]*a[2][0]*invdet;
          a1[1][0]=-a[0][1]*a[2][2]+a[0][2]*a[2][1]*invdet;
          a1[1][1]=a[0][0]*a[2][2]-a[0][2]*a[2][0]*invdet;
          a1[1][2]=-a[0][0]*a[2][1]+a[0][1]*a[2][0]*invdet;
          a1[2][0]=a[0][1]*a[1][2]-a[0][1]*a[1][1]*invdet;
          a1[2][1]=-a[0][0]*a[1][2]+a[0][2]*a[2][0]*invdet;
          a1[2][2]=a[0][0]*a[1][1]-a[0][1]*a[1][0]*invdet;
         }
        }
     //Ìàññèâ ôóíêöèé
        for (j=0; j<n; j++){
         f[0]=f1(x,y,z);
         f[1]=f2(x,y,z);
         f[2]=f3(x,y,z);
        }
     //Ìàññèâ ïåðåìåííûõ
        for (j=0; j<n; j++){
         K[0]=x;
         K[1]=y;
         K[2]=z;
        }
 
               do
                {
               for (i=0; i<n; i++){
                for (j=0; j<n; j++){
                 b[0]=a[0][0]*f[0]+a[0][1]*f[1]+a[0][2]*f[2];
                 b[1]=a[1][0]*f[0]+a[1][1]*f[1]+a[1][2]*f[2];
                 b[2]=a[2][0]*f[0]+a[2][1]*f[1]+a[2][2]*f[2];
                }
                }
                xk=K[0]-b[0];
                yk=K[1]-b[1];
                zk=K[2]-b[2];
 
                xn=xk-x;
                yn=yk-y;
                zn=zk-z;
                {
                 if (abs(xn)>abs(yn)){
                   if (abs(xn)>abs(zn))
                      max=xn;
                   else
                      max=zn;}
                 else {
                   if (abs(yn)>abs(zn))
                      max=yn;
                   else
                      max=zn;}
                   }
 
 
                 x=xk;
                 y=yk;
                 z=zk;}
                 while (abs(max)>=esp);
 
                 k++;
Edit5->Text=FloatToStr(xk);
Edit6->Text=FloatToStr(yk);
Edit7->Text=FloatToStr(zk);
Edit8->Text=FloatToStr(k);
}
}