import java.io.*;

/**
 * エムデン方程式を解くプログラム。
 */

public class Emdenequation
{
        static double n;

        public static void main()
        {
        double t,t0,x,v,dt,tmax;
        double k1[]=new double[2],k2[]=new double[2],k3[]=new double[2],k4[]=new double[2];
                String s;

                Emdenequation f = new Emdenequation();

                System.out.println("ポリトロープ指数の値は？");
                try
                {
                        InputStreamReader isr = new InputStreamReader(System.in);
                        BufferedReader br = new BufferedReader(isr);
                        s = br.readLine();
                        n = Double.parseDouble(s);
                        br.close();
                }
                catch (Exception e)
                {
                        System.out.println(e);
                }

                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;   //最大値

                try
                {
                        FileWriter fw = new FileWriter("output.data");  //gnuplot表示用データファイル

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

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

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

                                k4[0] = dt * f.f1(t + dt, x + k3[0], v + k3[1]);
                                k4[1] = dt * f.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;

                                fw.write(t + " " + x + " " + v + "\n");


                                if (x < 0) break;
                        }

                }
                catch(Exception e)
                {
                        System.out.print(e);
                }
        }
        double f1(double t, double x, double v)
        {
                return v;
        }

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

以下に上記のプログラムで求めた解を示す。横軸はξで縦軸はθである。ポリトロピック指数は0とした。