请各位学长帮忙修改个程序！解决思路

请各位学长帮忙修改个程序！
下面这个程序是用龙格库塔法求解微分方程的，哪位学长能帮忙把它改成用数组形式的，我不太会，想学习一下，我编的这个太繁杂可，谢谢！
/*龙格库塔法求解微分方程*/
#include "stdio.h"
#include "math.h"
#define A 0 /*定义初值和步长*/
#define B 0.002
#define C 1
#define D 1
#define E 1

double f(double,double,double,double); /* 定义函数*/
double g(double,double,double,double);
double e(double,double,double,double);

main()
{
  double x0=A, h=B,
  y01=C, u01=D, v01=E,
  j01,j02,j03,j04,k01,k02,k03,k04,l01,l02,l03,l04;
  int i,n=10000;
  FILE *cfPtr;
  if((cfPtr=fopen("luolunzi.dat","w"))==NULL) /*数据要写入的文件*/
  printf("File could not be opened\n");
  else
  { 
  for (i=1;i<=n;i++)
  { /*套用公式开始计算*/ 
  x0=x0+h;
  j01=f(x0,y01,u01,v01);
  k01=g(x0,y01,u01,v01);
  l01=e(x0,y01,u01,v01);
   
  j02=f(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);
  k02=g(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);
l02=e(x0+h/2.0,y01+h/2.0*j01,u01+h/2.0*k01,v01+h/2.0*l01);

  j03=f(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
  k03=g(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
  l03=e(x0+h/2.0,y01+h/2.0*j02,u01+h/2.0*k02,v01+h/2.0*l02);
   
  j04=f(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
  k04=g(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
  l04=e(x0+h,y01+h*j02,u01+h*k02,v01+h*l02);
   
  y01=y01+h/6.0*(j01+2.0*j02+2.0*j03+j04);
  u01=u01+h/6.0*(k01+2.0*k02+2.0*k03+k04);
  v01=v01+h/6.0*(l01+2.0*l02+2.0*l03+l04);
   
  printf("%15.7f\t%15.7f\t%15.7f\t\n",x0,y01,u01,v01);
  fprintf(cfPtr,"%15.7f\t%15.7f\t%15.7f\t\n",x0,y01,u01,v01);
  }
  fclose(cfPtr);
}
} 

double f(double s1,double q1,double u1,double r1) /*定义的微分方程*/  
{
  double z1;
  z1=-10*q1+(10*u1);
  return z1;
}
double g(double s2,double q2,double u2,double r2)
{
  double z2;
  z2=28*q2-u2-q2*r2;
  return z2;
}
double e(double s3,double q3,double u3,double r3)
{
double z3;
z3=(-8/3)*r3+q3*u3;
return z3;
}


------解决方案--------------------
楼主的意思是把每次计算出来的x0,y01,u01,v01用数组的形式保存下来,是吗?
------解决方案--------------------
整理修改如下：

C/C++ code

#include   <stdio.h> 
#include   <math.h> 

#define   A   0                                       /*定义初值和步长*/ 
#define   B   0.002 
#define   C   1 
#define   D   1 
#define   E   1 
#define   N   100

double   f1(double   s1,double   q1,double   u1,double   r1)         //f     
{ 
    return (-10*q1 + 10*u1); 
} 
double   f2(double   s2,double   q2,double   u2,double   r2)          //g
{ 
    return (28*q2 - u2 - q2*r2);
} 
double   f3(double   s3,double   q3,double   u3,double   r3)          //e
{ 
    return ((-8/3.0)*r3 + q3*u3);
} 

int main()
{
    double result[3] = {C, D, E};    //y,u,v
    double tmp[3][4];   //j,k,i
    //初始化
    double x = A, h = B;
    double (* p[3])(double, double, double, double) = {f1, f2, f3};
    
    int  i, j; 
    FILE *cfPtr; 
    
    if((cfPtr = fopen("luolunzi.dat","w")) == NULL)
        printf("File could not be opened\n"); 
    else 
    {   
        for(i = 1; i <= N; i++) 
          {         
            x += h;
            for(j = 0; j < 3; j++)
            {
                tmp[j][0] = p[j](x, result[0], result[1], result[2]);
            }    
            for(j = 0; j < 3; j++)
            {
                tmp[j][1] = p[j](x + h/2.0, result[0]+h/2.0*tmp[0][0], result[1]+h/2.0*tmp[1][0], result[2]+h/2.0*tmp[2][0]);
            }    
            for(j = 0; j < 3; j++)
            {
                tmp[j][2] = p[j](x + h/2.0, result[0]+h/2.0*tmp[0][1], result[1]+h/2.0*tmp[1][1], result[2]+h/2.0*tmp[2][1]);
            }
            for(j = 0; j < 3; j++)
            {
                tmp[j][3] = p[j](x + h, result[0]+h*tmp[0][1], result[1]+h*tmp[1][1], result[2]+h*tmp[2][1]);
            }
            printf("%15.7f\t", x);
            for(j = 0; j < 3; j++)
            {
                result[j] += h/6.0*(tmp[j][0]+2.0*tmp[j][1]+2.0*tmp[j][2]+tmp[j][3]);
                printf("%15.7f\t", result[j]);
            }
            fprintf(cfPtr,"%15.7f\t%15.7f\t%15.7f\t%15.7f\n", x, result[0], result[1], result[2]);
        }
        fclose(cfPtr);
    }
    return 0;
}