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

const int N=8;  //порядок точности
    double p[3][N+1];
    double t[3][N+1];
    double rad;
    double Ds;
    double Fi;
    double x=0;
    double y=0;
    
    double a=-sqrt(0.5);//окружность в пространстве
    double b=0;
    double c=sqrt(0.5);
    double x2=0;
    double y2=0;
    
    double dx=0;
    double dy=0;
    
    double kadr[2][7];
    double Pi=3.14159265359;
    double del=0; //ошибка
 
    p[0][0]=0;  // x0
    p[1][0]=rad;// y0
    p[2][0]=0;  // z0
    
 
/*
    double a=0;         //окружность в плоскости x-y
    double b=0;
    double c=1;
 
    p[0][0]=rad;// x0
    p[1][0]=0;  // y0
    p[2][0]=0;  // z0
*/
    t[0][0]=b*p[2][0]-c*p[1][0];  // f0=b*z0-c*y0
    t[1][0]=c*p[0][0]-a*p[2][0];  // g0=c*x0-a*z0
    t[2][0]=a*p[1][0]-b*p[0][0];  // h0=a*y0-b*x0
 
    // x0,y0,z0 - начальные значения координат;
    // a,b,c    - направляющие косинусы нормали
    //            к плоскости интерполяции.
 
    
 
    ofstream fout("sfera2.txt");
    cout.precision(w);
    cout.setf(ios::left,ios::adjustfield);
 
    cout << setw(v)<<"nс"<<setw(v)<<"x"<<setw(v)<<"y"
         <<setw(v)<<"z"<<setw(v)<<"del"<< endl;
 
    cout << setw(v)<<nс<<setw(v)<<p[0][0]<<setw(v)<<p[1][0]
         <<setw(v)<<p[2][0]<<setw(v)<<del<<endl;
 
    fout << setw(v)<<nс<<setw(v)<<p[0][0]<<setw(v)<<p[1][0]
             <<setw(v)<<p[2][0]<<setw(v)<<del<<endl;
 
//                      xn              yn                     xk               yk                  rad             fi
    kadr[0][0]=1, kadr[0][1]=100000, kadr[0][2]=0,          kadr[0][4]=0, kadr[0][5]=100000, kadr[0][7]=100000, kadr[0][8]=Pi/2; 
    kadr[1][0]=2, kadr[1][1]=0,      kadr[1][2]=-50000,     kadr[1][4]=0, kadr[1][5]=50000,  kadr[1][7]=50000,  kadr[1][8]=-Pi; 
 
    for (int k=0;k<2;k++)
    {
     rad=kadr[k][7];
     p[0]=kadr[k][1];//x0
     t[0]=kadr[k][2];//y0
 
     x1=p[0];
     y1=t[0];
     
 
     Ds=kadr[k][7]*rad/20;
     Fi=1/rad;
 
       for(int m=0;m<20;m++)
       {
        nс=nс+1;
        for (int i=1;i<N+1;i++)     //i-номер столбца
        {
            for (int j=1;j<N+1;j++) //j-номер строки
            {
              p[i][j]=Fi*t[i][j-1]/j;
              t[i][j]=-Fi*p[i][j-1]/j;
            };                      // кон.j
 
            for(int j=N; j>0; j=j-1)
            {
                p[i][j-1]=p[i][j]*Ds+p[i][j-1];
                t[i][j-1]=t[i][j]*Ds+t[i][j-1];
            };  
    
     };                         //кон. i
 
        del=sqrt(p[0][0]*p[0][0]+p[1][0]*p[1][0]+p[2][0]*p[2][0])-rad;
 
        x2=p[0],  y2=t[0];          //вычисление приращений и текущих 
        dx=x2-x1, dy=y2-y1;         //координат траектории обработки
        x1=x2,    y1=y2;
        x=x+dx,   y=y+dy;
 
        fout << setw(v)<<nс<<setw(v)<<p[0][0]<<setw(v)<<p[1][0]
             <<setw(v)<<p[2][0]<<setw(v)<<del<<endl;
 
        };                          //кон. m
        cout << setw(v)<<nс<<setw(v)<<p[0][0]<<setw(v)<<p[1][0]
             <<setw(v)<<p[2][0]<<setw(v)<<del<<endl;
 
    };                              //кон. к
    
    fout. close();
    return 0;
}