double** RKSystem4th( void fun(double,double[],double []), int N, double t0, double y00[],double tn, int M)
{
double** y=Make2DArray(M+1,N);
double k1[N],k2[N],k3[N],k4[N],temp[N];
double h=(tn-t0)/M;
int i=0,j=0;
y[0]=y00;
for(i=1;i<=2;i++)
{
fun(t0,y00,k1);
printf("///////y00=%lf,%lf\n",y00[0],y00[1]);
arrayadd1(y00,k1,temp,h/2,N);
printf("///////y00=%lf,%lf\n",y00[0],y00[1]);
fun(t0+h/2,temp,k2);
arrayadd1(y00,k2,temp,h/2,N);
fun(t0+h/2,temp,k3);
arrayadd1(y00,k3,temp,h,N);
fun(t0+h,temp,k4);
arrayadd1(k1,k2,temp,2,N);
arrayadd1(temp,k3,temp,2,N);
arrayadd1(temp,k4,temp,1,N);
arrayadd1(y00,temp,temp,h/6,N);
y[i]=temp;
y00=temp;printf("y00=%lf,%lf\n\n",y00[0],y00[1]);
t0=t0+h;
}
return y;
}
void arrayadd1(const double a[],double b[], double c[],double h,int n)
{
int i;
for(i=0;i<=n-1;i++)
{
c[i]=a[i]+h*b[i];
}
}