// ================================================================================
// 以下、4次ルンゲ=クッタ(RK4)法で解曲線の座標値を計算
// ================================================================================

for( int i=0; i<LOOPS; i++ ){

	// グラフプロット周期ごとに座標値を書き出す
	if( i % OUTPUT_INTERVAL == 0 ){
		println(x, y, z);
	}

	double k1_x = dxdt(x, y, z);
	double k1_y = dydt(x, y, z);
	double k1_z = dzdt(x, y, z);

	double k2_x = dxdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
	double k2_y = dydt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);
	double k2_z = dzdt(x + k1_x*DT/2.0, y + k1_y*DT/2.0, z + k1_z*DT/2.0);

	double k3_x = dxdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
	double k3_y = dydt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);
	double k3_z = dzdt(x + k2_x*DT/2.0, y + k2_y*DT/2.0, z + k2_z*DT/2.0);

	double k4_x = dxdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
	double k4_y = dydt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);
	double k4_z = dzdt(x + k3_x*DT, y + k3_y*DT, z + k3_z*DT);

	x += (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x) * DT / 6.0;
	y += (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y) * DT / 6.0;
	z += (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z) * DT / 6.0;
}

// 最後の座標値を書き出しておく
println(x, y, z);