// ==================================================================
// RungeKutta.sci: 4次のルンゲ・クッタ法(6.42)で微分方程式を解く関数
// ------------------------------------------------------------------
//  ＜引数＞     f: 微分方程式の関数 f(x,y)
//             a,b: 区間 [a,b]
//              y0: 初期条件
//               N: 区間の分割数
// ＜返り値＞  x,y: 微分方程式の数値解 (x:ベクトル, y:ベクトル)
// ==================================================================
function [x,y] = RungeKutta(f, a, b, y0, N)
h = (b - a)/N;                                    // h：区間の刻み幅
x = a:h:b;                                        // x軸の刻み(縦ベクトル)
y = y0 ;                                          // 初期条件y=y0
for i = 1 : N                                     // N点で解を求める
   k1 = h * f( x(i),       y(i)        );         // 式(6.42) の k1
   k2 = h * f( x(i) + h/2, y(i) + k1/2 );         // 式(6.42) の k2
   k3 = h * f( x(i) + h/2, y(i) + k2/2 );         // 式(6.42) の k3
   k4 = h * f( x(i) + h,   y(i) + k3   );         // 式(6.42) の k4
   y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4) / 6;   // 式(6.42) の y(i+1)
end
endfunction