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

double n;       //ポリトロピック指数

double f1(double t,double x,double v);
double f2(double t,double x,double v);

int main()
{
        double t,t0,x,v,dt,tmax;
        double k1[2],k2[2],k3[2],k4[2];

                t0=0.0001;      //初期値
        x=1.0-1.0/6.0*t0*t0+n/120.0*t0*t0*t0*t0;        //初期値
        v=-1.0/3.0*t0+n/30.0*t0*t0*t0;  //初期値
        dt=0.001;       //刻み幅
        tmax=10000.0;   //最大値

        FILE *output,*outputt,*outputx,*outputv;
        output=fopen("output.data","w");        //gnuplot表示用データファイル
                outputt=fopen("outputt.data","w");      //excel表示用データファイル
                outputx=fopen("outputx.data","w");
                outputv=fopen("outputv.data","w");

                printf("ポリトロープ指数の値は？");
                scanf("%lf",&n);

        for(t=t0;t<tmax;t+=dt) {
                k1[0]=dt*f1(t,x,v);
                k1[1]=dt*f2(t,x,v);

                k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);
                k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);

                k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);
                k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);

                k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]);
                k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]);

                x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0;
                v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0;

                fprintf(output,"%f %f %f\n",t,x,v);
                                fprintf(outputt,"%f\n",t);
                                fprintf(outputx,"%f\n",x);
                                fprintf(outputv,"%f\n",v);

                                if(x<0) break;
        }

        fclose(output);

        return 0;
}

double f1(double t,double x,double v)
{
        return v;
}

double f2(double t,double x,double v)
{
        return -pow(x,n)-2.0/t*v;
}

以下に上記のプログラムで求めた解を示す。グラフにはExcel(エクセル)を用いた。