大神们看下这个程序怎么改

佩刀

#include<iostream>
#include<iomanip>
using namespace std;
void RK4(double(*f)(double t,double x,double y,double z),double(*g)(double t,double x,double y,double z),double(*k)(double t,double x,double y,double z),double initial[4],double resu[4],double h)
{
double f1,f2,f3,f4,g1,g2,g3,g4,k1,k2,k3,k4,t0,x0,y0,z0,x1,y1,z1;
t0=initial[0];x0=initial[1];y0=initial[2];z0=initial[3];

f1=f(t0,x0,y0,z0);                                   g1=g(t0,x0,y0,z0);                                     k1=k(t0,x0,y0,z0);
f2=f((t0+h/2),(x0+h*f1/2),(y0+h*g1/2),(z0+h*k1/2));  g2=g((t0+h/2),(x0+h*f1/2),(y0+h*g1/2),(z0+h*k1/2));    k2=k((t0+h/2),(x0+h*f1/2),(z0+h*k1/2));
f3=f((t0+h/2),(x0+h*f2/2),(y0+h*g2/2),(z0+h*k2/2));  g3=g((t0+h/2),(x0+h*f2/2),(y0+h*g2/2),(z0+h*k2/2));    k3=k((t0+h/2),(x0+h*f2/2),(z0+h*k2/2));
f4=f((t0+h),(x0+h*f3),(y0+h*g3),(z0+h*k3));          g4=((t0+h),(x0+h*f3),(y0+h*g3),(z0+h*k3));             k4=k((t0+h),(x0+h*f3),(z0+h*k3));

x1=x0+h*(f1+2*f2+2*f3+f4)/6;      
y1=y0+h*(g1+2*g2+2*g3+g4)/6;      
z1=z0+h*(k1+2*k2+2*k3+k4)/6;
resu[0]=t0+h;resu[1]=x1;resu[2]=y1;resu[3]=z1;
}
int main()
{
double f(double t,double x,double y,double z);
double g(double t,double x,double y,double z);
double k(double t,double x,double y,double z);
double initial[4],resu[4];
double t,step;
double a,b,H;
int i;
cout<<"输入所求微分方程组的初值t0,x0,y0,z0:";
cin>>initial[0]>>initial[1]>>initial[2]>>initial[3];
cout<<"输入所求微分方程组的微分区间[a,b]:";
cin>>a>>b;
cout<<"输入所求微分方程组所分解子区间的个数step:";
cin>>step;
cout<<setiosflags(ios::right)<<setiosflags<<(ios::fixed)<<setiosflags(10);
H=(b-a)/step;
cout<<initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<setw(18)<<initial[3]<<endl;
for(i=0;i<step;i++)
{
RK4(f,g,k,initial,resu,H);
cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<setw(20)<<resu[3]<<endl;
initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];initial[3]=resu[3];
}
return(0);
}

double f(double t,double x,double y,double z)
{
double dx;
dx=-10*x+10*y;
return(dx);
}
double g(double t,double x,double y,double z)
{
double dy;
dy=30*x-y-x*z;
return(dy);
}
double k(double t,double x,double y,double z)
{
double dz;
dz=-(8/3)*z+x*y;
return(dz);
}