説明できるけどテキストベースじゃ面倒なので省略 
続き 
ルンゲクッタ−マースン法 
ルンゲクッタ法は与えられた関数が素直な場合は問題ないが 
解に不連続点を持つ物や関数の形が複雑な物は計算できない 
そこでルンゲクッタ法の基本式を次のように変形する 
k0=h*(f(xi,yi)/3 
k1=h*f(xi+h/3,yi+k0)/3 
k2=h*f(xi+h/3,yi+k0/2+k1/2)/3 
k3=h*f(xi+h/2,yi+3/8*k0+9/8*k2)/3 
k4=h*f(xi+h,yi+3/2*k0-9/2*k2+6*k3)/3 
yi+1=yi+(k0+4*k3+k4)/2+O(h^5) 
ここで許容誤差をεとした場合絶対誤差ε'は以下の式で計算できる 
ε'=(k0-9/2*k2+4*k3-1/2*k4)/5 
εとε'との関係は計算刻み幅hと関連し 
|ε'|＜εならばyi+1は計算できたことになるが特に 
|ε'|＜ε/32ならば刻み幅を２倍にできる 
|ε'|≧ε/32ならば同じ幅で計算 
|ε'|≧εならば今の計算幅が不適当だったので刻み幅をh=h/2 
にして再度ε'を計算し|ε'|＜εの条件が成立するまで繰り返す 

public double[] RungeKuttaGill(double Y0,double XST,double DX,int N)
   { double[]Y = new double[N]; double X; int i; double S2=Math.sqrt(2);
     double A1 = 1 / 4 - 1 / (2*S2);  double A2 = 1/2 - 1 /(2* S2);
     double B1 = 1 / S2            ;  double B2 = 1 + 1 / S2;
     double C1 =  2 - 2 / S2       ;  double C2 = 2 + 2 / S2;
  
     for(Y[0] = Y0, X = XST, i = 1; i < N; i++, X += DX)
     {
       double K1 = g(X, Y[i - 1]);
       double K2 = g(X + DX / 2, Y[i - 1] + DX * K1 / 2                );
       double K3 = g(X + DX / 2, Y[i - 1] - A1 * DX * K1 + A2 * DX * K2);
       double K4 = g(X + DX    , Y[i - 1] - B1 * DX * K2 + B2 * DX * K3);
       Y[i] = Y[i - 1] + (K1 + C1 * K2 + C2 * K3 + K4) * DX / 6;
     }
     return Y;
   }
