#include "stdio.h"
#include "stdlib.h"

void RKT(t,y,n,h,k,z)
int n;				/*微分方程組中方程的個數(shù)，也是未知函數(shù)的個數(shù)*/
int k;				/*積分的步數(shù)(包括起始點這一步)*/
double t;			/*積分的起始點t0*/
double h;			/*積分的步長*/
double y[];			/*存放n個未知函數(shù)在起始點t處的函數(shù)值,返回時,其初值在二維數(shù)組z的第零列中*/
double z[];			/*二維數(shù)組,體積為n x k.返回k個積分點上的n個未知函數(shù)值*/
{
	extern void Func();				/*聲明要求解的微分方程組*/
    int i,j,l;
    double a[4],*b,*d;
    b=malloc(n*sizeof(double));		/*分配存儲空間*/
	if(b == NULL)
	{
		printf("內(nèi)存分配失敗\n");
		exit(1);
	}
    d=malloc(n*sizeof(double));		/*分配存儲空間*/
	if(d == NULL)
	{
		printf("內(nèi)存分配失敗\n");
		exit(1);
	}
	/*后面應(yīng)用RK4公式中用到的系數(shù)*/
    a[0]=h/2.0;						
	a[1]=h/2.0;
    a[2]=h; 
	a[3]=h;
    for(i=0; i<=n-1; i++) 
		z[i*k]=y[i];				/*將初值賦給數(shù)組z的相應(yīng)位置*/
    for(l=1; l<=k-1; l++)
    {
		Func(y,d);
        for (i=0; i<=n-1; i++)
			b[i]=y[i];
        for (j=0; j<=2; j++)
        {
			for (i=0; i<=n-1; i++)
            {
				y[i]=z[i*k+l-1]+a[j]*d[i];
                b[i]=b[i]+a[j+1]*d[i]/3.0;
            }
            Func(y,d);
        }
        for(i=0; i<=n-1; i++)
          y[i]=b[i]+h*d[i]/6.0;
        for(i=0; i<=n-1; i++)
          z[i*k+l]=y[i];
        t=t+h;
	}
    free(b);			/*釋放存儲空間*/
	free(d);			/*釋放存儲空間*/
    return;
}
main()
{
	int i,j;
    double t,h,y[3],z[3][11];
    y[0]=-1.0; 
	y[1]=0.0; 
	y[2]=1.0;
    t=0.0; 
	h=0.01;
    RKT(t,y,3,h,11,z);
    printf("\n");
    for (i=0; i<=10; i++)			/*打印輸出結(jié)果*/
    {
		t=i*h;
        printf("t=%5.2f\t   ",t);
        for (j=0; j<=2; j++)
          printf("y(%d)=%e  ",j,z[j][i]);
        printf("\n");
	}
}

void Func(y,d)
double y[],d[];
{
	d[0]=y[1];		/*y0'=y1*/
	d[1]=-y[0];		/*y1'=y0*/
	d[2]=-y[2];		/*y2'=y2*/
	return;
}