这个程序的A[k+1]怎么会是9562000000000000000000.000000，哪里错了呢？
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
#include <iostream.h>
void main()
{
double A[100],f[1000],p[1000],q[1000],x[1000],y[1000],z[1000],l[1000],m[1000],n[1000];
int i,j,k;
double F1,F2,F3,F4,P1,P2,P3,P4,Q1,Q2,Q3,Q4;
double X1,X2,X3,X4,Y1,Y2,Y3,Y4,Z1,Z2,Z3,Z4;
double L1,L2,L3,L4,M1,M2,M3,M4,N1,N2,N3,N4;
double h=0.05;
/*printf("请输入第一个A的值\n");
cin>>A[0];*/
A[0]=0.2;
for(k=0;k<100;k++)
{
   
    for(j=0;j<1000;j++)
    {
        f[0]=0,p[0]=0,q[0]=A[k];

        F1=h*p[j];
        P1=h*q[j];
        Q1=h*(p[j]*p[j]-f[j]*q[j]-1);

        F2=h*(p[j]+0.5*P1);
        P2=h*(q[j]+0.5*Q1);
        Q2=h*((p[j]+0.5*P1)*(p[j]+0.5*P1)-(f[j]+0.5*F1)*(q[j]+0.5*Q1)-1);

        F3=h*(p[j]+0.5*P2);
        P3=h*(q[j]+0.5*Q2);
        Q3=h*((p[j]+0.5*P2)*(p[j]+0.5*P2)-(f[j]+0.5*F2)*(q[j]+0.5*Q2)-1);

        F4=h*(p[j]+P3);
        P4=h*(q[j]+Q3);
        Q4=h*((p[j]+P3)*(p[j]+P3)-(f[j]+F3)*(q[j]+Q3)-1);

        f[j+1]=f[j]+(1/6.0)*(F1+2*F2+2*F3+F4);
        p[j+1]=p[j]+(1/6.0)*(P1+2*P2+2*P3+P4);
        q[j+1]=p[j]+(1/6.0)*(Q1+2*Q2+2*Q3+Q4);

        if(p[j+1]-p[j]<1e-10) break;
        
    }

    if((p[j+1]-1)<1e-6) break;
   


    x[0]=0,y[0]=0,z[0]=A[k],l[0]=0,m[0]=0,n[0]=0;
    for(i=0;i<1000;i++)
    {
        X1=h*y[i];
        Y1=h*z[i];
        Z1=h*(y[i]*y[i]-x[i]*z[i]-1);
        L1=h*m[i];
        M1=h*n[i];
        N1=h*(-z[i]*l[i]+2*y[i]*m[i]-x[i]*n[i]);

        X2=h*(y[i]+0.5*Y1);
        Y2=h*(z[i]+0.5*Z1);
        Z2=h*((y[i]+0.5*Y1)*(y[i]+0.5*Y1)-(x[i]+0.5*X1)*(z[i]+0.5*Z1)-1);
        L2=h*(m[i]+0.5*M1);
        M2=h*(n[i]+0.5*N1);
        N2=h*(-(z[i]+0.5*Z1)*(l[i]+0.5*L1)+2*(y[i]+0.5*Y1)*(m[i]+0.5*M1)-(x[i]+0.5*X1)*(n[i]+0.5*N1));

        X3=h*(y[i]+0.5*Y2);
        Y3=h*(z[i]+0.5*Z2);
        Z3=h*((y[i]+0.5*Y2)*(y[i]+0.5*Y2)-(x[i]+0.5*X2)*(z[i]+0.5*Z2)-1);
        L3=h*(m[i]+0.5*M2);
        M3=h*(n[i]+0.5*N2);
        N3=h*(-(z[i]+0.5*Z2)*(l[i]+0.5*L2)+2*(y[i]+0.5*Y2)*(m[i]+0.5*M2)-(x[i]+0.5*X2)*(n[i]+0.5*N2));

        X4=h*(y[i]+Y3);
        Y4=h*(z[i]+Z3);
        Z4=h*((y[i]+Y3)*(y[i]+Y3)-(x[i]+X3)*(z[i]+Z3)-1);
        L4=h*(m[i]+M3);
        M4=h*(n[i]+N3);
        N4=h*(-(z[i]+Z3)*(l[i]+L3)+2*(y[i]+Y3)*(m[i]+M3)-(x[i]+X3)*(n[i]+N3));

        x[i+1]=x[i]+(1/6.0)*(X1+2*X2+2*X3+X4);
        y[i+1]=y[i]+(1/6.0)*(Y1+2*Y2+2*Y3+Y4);
        z[i+1]=z[i]+(1/6.0)*(Z1+2*Z2+2*Z3+Z4);
        l[i+1]=l[i]+(1/6.0)*(L1+2*L2+2*L3+L4);
        m[i+1]=m[i]+(1/6.0)*(M1+2*M2+2*M3+M4);
        n[i+1]=n[i]+(1/6.0)*(N1+2*N2+2*N3+N4);

        if((y[i+1]-y[i])<1e-6&&(m[i+1]-m[i])<1e-6) break;


    }
        A[k+1]=A[k]-(y[i+1]-1)/m[i+1];

}

printf("打靶修正后的A=%8.6f\n",A[k+1]);

}