微分方程式（ルンゲ・クッタ）

　　プログラムは，以下の微分方程式をルンゲ・クッタ法によって，0 から 1 秒まで解いた例です．

　　d²y/dt² + 3dy/dt + 2y = 1　　初期条件はすべて０

　　（x[0] = y, x[1] = dy/dt)

/**********************************/
/* 微分方程式（ルンゲ・クッタ法） */
/*      例：d2x/dt+3dx/dt+2x=1    */
/*      coded by Y.Suganuma       */
/**********************************/
#include <stdio.h>

void snx(double, double *, double *);
double rkg(double, double, double *, double *, double **, int,
           void (*)(double, double *, double *));
int main()
{
	double time, h, *x, *dx, **g;
	int i1, n = 2;
/*
          初期設定
*/
	x  = new double [n];
	dx = new double [n];
	g  = new double * [4];
	for (i1 = 0; i1 < 4; i1++)
		g[i1] = new double [n];
	time = 0.0;
	h    = 0.01;
	x[0] = 0.0;
	x[1] = 0.0;
/*
          計算と出力
*/
	for (i1 = 0; i1 < 101; i1++) {
		printf("time = %f, x = %f\n", time, x[0]);
		time = rkg(time, h, x, dx, g, n, snx);
	}

	return 0;
}

/****************/
/* 微係数の計算 */
/****************/
void snx(double time, double *x, double *dx)
{
	dx[0] = x[1];
	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
}

/*******************************************/
/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
/*      time : 現在の時間                  */
/*      h : 時間刻み幅                     */
/*      x : 現在の状態                     */
/*      dx : 微係数（f(t,x):snxで計算する）*/
/*      g : 作業域（g[4][n]）              */
/*      n : 微分方程式の次数               */
/*      sub : 微係数を計算する関数の名前   */
/*      return : time+h                    */
/*******************************************/
double rkg(double time, double h, double *x, double *dx, double **g,
           int n, void (*sub)(double, double *, double *))
{
	int i1;
	double h2;

	h2 = 0.5 * h;

	(*sub)(time, x, dx);
	for (i1 = 0; i1 < n; i1++)
		g[0][i1] = h * dx[i1];

	time += h2;
	for (i1 = 0; i1 < n; i1++)
		g[1][i1] = x[i1] + 0.5 * g[0][i1];
	(*sub)(time, g[1], dx);
	for (i1 = 0; i1 < n; i1++)
		g[1][i1] = h * dx[i1];

	for (i1 = 0; i1 < n; i1++)
		g[2][i1] = x[i1] + 0.5 * g[1][i1];
	(*sub)(time, g[2], dx);
	for (i1 = 0; i1 < n; i1++)
		g[2][i1] = h * dx[i1];

	time += h2;
	for (i1 = 0; i1 < n; i1++)
		g[3][i1] = x[i1] + g[2][i1];
	(*sub)(time, g[3], dx);
	for (i1 = 0; i1 < n; i1++)
		g[3][i1] = h * dx[i1];

	for (i1 = 0; i1 < n; i1++)
		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;

	return time;
}

Java

/**********************************/
/* 微分方程式（ルンゲ・クッタ法） */
/*      例：d2x/dt+3dx/dt+2x=1    */
/*      coded by Y.Suganuma       */
/**********************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double time, h, x[], dx[], g[][];
		int i1, n = 2;
	/*
	          初期設定
	*/
		x    = new double [n];
		dx   = new double [n];
		g    = new double [4][n];

		time = 0.0;
		h    = 0.01;
		x[0] = 0.0;
		x[1] = 0.0;
	/*
	          計算と出力
	*/
		Kansu kn = new Kansu();

		for (i1 = 0; i1 < 101; i1++) {
			System.out.println("time = " + time + " x = " + x[0]);
			time = kn.rkg(time, h, x, dx, g, n);
		}
	}
}

/****************/
/* 微係数の計算 */
/****************/
class Kansu extends RKG
{
	void snx(double time, double x[], double dx[])
	{
		dx[0] = x[1];
		dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
	}
}

abstract class RKG
{
	/*************************************************/
	/* ルンゲ・クッタ法  dx/dt=f(t,x)                */
	/*      time : 現在の時間                        */
	/*      h : 時間刻み幅                           */
	/*      x : 現在の状態                           */
	/*      dx : 微係数（f(t,x):snxで計算する）      */
	/*      g : 作業域（g[4][n]）                    */
	/*      n : 微分方程式の次数                     */
	/*      return : time+h                          */
	/*************************************************/

	abstract void snx(double time, double x[], double dx[]);

	double rkg(double time, double h, double x[], double dx[], double g[][], int n)
	{
		int i1;
		double h2;

		h2 = 0.5 * h;

		snx(time, x, dx);
		for (i1 = 0; i1 < n; i1++)
			g[0][i1] = h * dx[i1];

		time += h2;
		for (i1 = 0; i1 < n; i1++)
			g[1][i1] = x[i1] + 0.5 * g[0][i1];
		snx(time, g[1], dx);
		for (i1 = 0; i1 < n; i1++)
			g[1][i1] = h * dx[i1];

		for (i1 = 0; i1 < n; i1++)
			g[2][i1] = x[i1] + 0.5 * g[1][i1];
		snx(time, g[2], dx);
		for (i1 = 0; i1 < n; i1++)
			g[2][i1] = h * dx[i1];

		time += h2;
		for (i1 = 0; i1 < n; i1++)
			g[3][i1] = x[i1] + g[2][i1];
		snx(time, g[3], dx);
		for (i1 = 0; i1 < n; i1++)
			g[3][i1] = h * dx[i1];

		for (i1 = 0; i1 < n; i1++)
			x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;

		return time;
	}
}

JavaScript

　　ここをクリックすると，JavaScript の仕様に適合した形で微分方程式を入力することによって，任意の微分方程式の解を画面上で求めることができます．

<!DOCTYPE HTML>

<HTML>

<HEAD>

	<TITLE>微分方程式（ルンゲ・クッタ）</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		eq = new Array();

		function main()
		{
					// データの設定
			let n  = parseInt(document.getElementById("order").value);
			let t1 = parseFloat(document.getElementById("t1").value);
			let t2 = parseFloat(document.getElementById("t2").value);
			let h  = parseFloat(document.getElementById("haba").value);
			let x  = new Array();
			x = document.getElementById("ini").value.split(/ {1,}/)
			for (let i1 = 0; i1 < n; i1++)
				x[i1] = parseFloat(x[i1]);
			eq = document.getElementById("func").value.split(/\n{1,}/)
			for (let i1 = 0; i1 < n; i1++)
				eq[i1] = eq[i1] + ";";
			dx = new Array();
			g  = new Array();
			for (let i1 = 0; i1 < 4; i1++)
				g[i1] = new Array();
					// 実行，結果
			let time = t1;
			let c = 100000;
			let str  = "";
			while (time < t2+0.1*h) {
				let s1 = Math.round(time * c) / c;
				str = str + "time = " + s1 + " x =";
				for (let i1 = 0; i1 < n; i1++) {
					let s2 = Math.round(x[i1] * c) / c;
					str = str + " " + s2;
				}
				str += "\n";
				time = rkg(time, h, x, dx, g, n, snx);
			}
			document.getElementById("res").value = str;
		}

		/****************/
		/* 関数値の計算 */
		/****************/
		function snx(n, x, dx)
		{
			for (let i1 = 0; i1 < n; i1++)
				dx[i1] = eval(eq[i1]);
		}

		/********************************************/
		/* ルンゲ・クッタ法  dx/dt=f(t,x)           */
		/*      time : 現在の時間                   */
		/*      h : 時間刻み幅                      */
		/*      x : 現在の状態                      */
		/*      dx : 微係数（f(t,x):snxで計算する） */
		/*      g : 作業域（g[4][n]）               */
		/*      n : 微分方程式の次数                */
		/*      fn : 関数値を計算する関数　 　　　　*/
		/*      return : time+h                     */
		/********************************************/
		function rkg(time, h, x, dx, g, n, fn)
		{
			let h2 = 0.5 * h;

			fn(n, x, dx);
			for (let i1 = 0; i1 < n; i1++)
				g[0][i1] = h * dx[i1];

			time += h2;
			for (let i1 = 0; i1 < n; i1++)
				g[1][i1] = x[i1] + 0.5 * g[0][i1];
			fn(n, g[1], dx);
			for (let i1 = 0; i1 < n; i1++)
				g[1][i1] = h * dx[i1];

			for (let i1 = 0; i1 < n; i1++)
				g[2][i1] = x[i1] + 0.5 * g[1][i1];
			fn(n, g[2], dx);
			for (let i1 = 0; i1 < n; i1++)
				g[2][i1] = h * dx[i1];

			time += h2;
			for (let i1 = 0; i1 < n; i1++)
				g[3][i1] = x[i1] + g[2][i1];
			fn(n, g[3], dx);
			for (let i1 = 0; i1 < n; i1++)
				g[3][i1] = h * dx[i1];

			for (let i1 = 0; i1 < n; i1++)
				x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;

			return time;
		}
	</SCRIPT>

</HEAD>

<BODY STYLE="font-size: 130%; background-color: #eeffee;">

	<H2 STYLE="text-align:center"><B>微分方程式（ルンゲ・クッタ）</B></H2>

	<DL>
		<DT>　　テキストフィールド及びテキストエリアには，例として，微分方程式，
		<P STYLE="text-align:center"><IMG SRC="rkg.gif"></P>
		<DT>を解く場合に対する値が設定されています（他の問題を実行する場合は，それらを適切に修正してください）．ルンゲ・クッタ法を使用するためには，1 次の連立微分方程式に変換する必要があるため，
		<P STYLE="text-align:center"><IMG SRC="rkg1.gif"></P>
		<DT>とおき，以下に示すような 1 次の連立微分方程式に変換しています．
		<P STYLE="text-align:center"><IMG SRC="rkg2.gif"></P>
		<DT>テキストエリアに設定してあるのは，これらの式の右辺です．ただし，変数 x<SUB>1</SUB> や x<SUB>2</SUB> の代わりに，配列の要素 x[0] と x[1] を使用しています．なお，式は，JavaScript の仕様に適合した形式で記述してあることに注意してください．
	</DL>

	<DIV STYLE="text-align:center">

		次数(n)：<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="2">　
		範囲：<INPUT ID="t1" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.0">，
				<INPUT ID="t2" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="1.0">　
		幅：<INPUT ID="haba" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.01">　
		初期値(n個)：<INPUT ID="ini" STYLE="font-size: 100%" TYPE="text" SIZE="10" VALUE="0.0 0.0"><BR><BR>
		微分方程式　　<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR>
		dx = <TEXTAREA ID="func" STYLE="font-size: 110%" COLS="70" ROWS="10">
x[1]
-2.0 * x[0] - 3.0 * x[1] + 1.0
		</TEXTAREA><BR><BR>
		結果<BR>
		<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="50" ROWS="10"></TEXTAREA>
	</DIV>

</BODY>

</HTML>

PHP

<?php

/*
          初期設定
*/
	$n  = 2;
	$x  = array($n);
	$dx = array($n);
	$g  = array(4);
	for ($i1 = 0; $i1 < 4; $i1++)
		$g[$i1] = array($n);
	$time = 0.0;
	$h    = 0.01;
	$x[0] = 0.0;
	$x[1] = 0.0;
/*
          計算と出力
*/
	for ($i1 = 0; $i1 < 101; $i1++) {
		printf("time = %f, x = %f\n", $time, $x[0]);
		$time = rkg($time, $h, $x, $dx, $g, $n, "snx");
	}

/****************/
/* 微係数の計算 */
/****************/
function snx($time, $x, &$dx)
{
	$dx[0] = $x[1];
	$dx[1] = -2.0 * $x[0] - 3.0 * $x[1] + 1.0;
}

/*******************************************/
/* ルンゲ・クッタ法  dx/dt=f(t,x)          */
/*      time : 現在の時間                  */
/*      h : 時間刻み幅                     */
/*      x : 現在の状態                     */
/*      dx : 微係数（f(t,x):snxで計算する）*/
/*      g : 作業域（g[4][n]）              */
/*      n : 微分方程式の次数               */
/*      sub : 微係数を計算する関数の名前   */
/*      return : time+h                    */
/*******************************************/
function rkg($time, $h, &$x, $dx, $g, $n, $sub)
{
	$h2 = 0.5 * $h;

	$sub($time, $x, $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[0][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $x[$i1] + 0.5 * $g[0][$i1];
	$sub($time, $g[1], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[1][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $x[$i1] + 0.5 * $g[1][$i1];
	$sub($time, $g[2], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[2][$i1] = $h * $dx[$i1];

	$time += $h2;
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $x[$i1] + $g[2][$i1];
	$sub($time, $g[3], $dx);
	for ($i1 = 0; $i1 < $n; $i1++)
		$g[3][$i1] = $h * $dx[$i1];

	for ($i1 = 0; $i1 < $n; $i1++)
		$x[$i1] = $x[$i1] + ($g[0][$i1] + 2.0 * $g[1][$i1] + 2.0 * $g[2][$i1] + $g[3][$i1]) / 6.0;

	return $time;
}

?>

Ruby

#*********************************/
# 微分方程式（ルンゲ・クッタ法） */
#      例：d2x/dt+3dx/dt+2x=1    */
#      coded by Y.Suganuma       */
#*********************************/

#***************/
# 微係数の計算 */
#***************/
snx = Proc.new { |time, x, dx|
	dx[0] = x[1]
	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0
}

#******************************************/
# ルンゲ・クッタ法  dx/dt=f(t,x)          */
#      time : 現在の時間                  */
#      h : 時間刻み幅                     */
#      x : 現在の状態                     */
#      dx : 微係数（f(t,x):snxで計算する）*/
#      g : 作業域（g[4][n]）              */
#      n : 微分方程式の次数               */
#      snx : 微係数を計算する関数の名前   */
#      return : time+h                    */
#******************************************/
def rkg(time, h, x, dx, g, n, &sub)

	h2 = 0.5 * h

	sub.call(time, x, dx)
	for i1 in 0 ... n
		g[0][i1] = h * dx[i1]
	end

	time += h2
	for i1 in 0 ... n
		g[1][i1] = x[i1] + 0.5 * g[0][i1]
	end
	sub.call(time, g[1], dx)
	for i1 in 0 ... n
		g[1][i1] = h * dx[i1]
	end

	for i1 in 0 ... n
		g[2][i1] = x[i1] + 0.5 * g[1][i1]
	end
	sub.call(time, g[2], dx)
	for i1 in 0 ... n
		g[2][i1] = h * dx[i1]
	end

	time += h2
	for i1 in 0 ... n
		g[3][i1] = x[i1] + g[2][i1]
	end
	sub.call(time, g[3], dx)
	for i1 in 0 ... n
		g[3][i1] = h * dx[i1]
	end

	for i1 in 0 ... n
		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0
	end

	return time
end

		# 初期設定
n  = 2
x  = Array.new(n)
dx = Array.new(n)
g  = Array.new(4)
for i1 in 0 ... 4
	g[i1] = Array.new(n)
end
time = 0.0
h    = 0.01
x[0] = 0.0
x[1] = 0.0
		# 計算と出力
for i1 in 0 ... 101
	printf("time = %f, x = %f\n", time, x[0])
	time = rkg(time, h, x, dx, g, n, &snx)
end

Python

# -*- coding: UTF-8 -*-
import numpy as np
from math import *

############################################
# ルンゲ・クッタ法  dx/dt=f(t,x)
#      time : 現在の時間
#      h : 時間刻み幅
#      x : 現在の状態
#      dx : 微係数（f(t,x):subで計算する）
#      g : 作業域（g[4][n]）
#      n : 微分方程式の次数
#      sub : 微係数を計算する関数の名前
#      return : time+h
#      coded by Y.Suganuma
############################################

def rkg(time, h, x, dx, g, n, sub) :

	h2 = 0.5 * h
	sub(time, x, dx)
	for i1 in range(0, n) :
		g[0][i1] = h * dx[i1]

	time += h2
	for i1 in range(0, n) :
		g[1][i1] = x[i1] + 0.5 * g[0][i1]
	sub(time, g[1], dx)
	for i1 in range(0, n) :
		g[1][i1] = h * dx[i1]

	for i1 in range(0, n) :
		g[2][i1] = x[i1] + 0.5 * g[1][i1]
	sub(time, g[2], dx)
	for i1 in range(0, n) :
		g[2][i1] = h * dx[i1]

	time += h2
	for i1 in range(0, n) :
		g[3][i1] = x[i1] + g[2][i1]
	sub(time, g[3], dx)
	for i1 in range(0, n) :
		g[3][i1] = h * dx[i1]

	for i1 in range(0, n) :
		x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0

	return time

############################################
# 微分方程式（ルンゲ・クッタ法）
#      例：d2x/dt+3dx/dt+2x=1
#      coded by Y.Suganuma
############################################

			# 微係数の計算 */
def snx(time, x, dx) :
	dx[0] = x[1]
	dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0
			# 初期設定
n    = 2
x    = np.array([0.0, 0.0], np.float)
dx   = np.empty(n, np.float)
g    = np.empty((4, n), np.float)
time = 0.0
h    = 0.01
			# 計算と出力
for i1 in range(0, 101) :
	print("time = " + str(time) + ", x = " + str(x[0]))
	time = rkg(time, h, x, dx, g, n, snx)

/**********************************/
/* 微分方程式（ルンゲ・クッタ法） */
/*      例：d2x/dt+3dx/dt+2x=1    */
/*      coded by Y.Suganuma       */
/**********************************/
using System;

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	public Test1()
	{
	/*
	          初期設定
	*/
		int n        = 2;
		double time  = 0.0;
		double h     = 0.01;
		double[] x   = {0.0, 0.0};
		double[] dx  = new double [n];
		double[][] g = new double [4][];
		g[0] = new double[2];
		g[1] = new double[2];
		g[2] = new double[2];
		g[3] = new double[2];
	/*
	          計算と出力
	*/
		for (int i1 = 0; i1 < 101; i1++) {
			Console.WriteLine("time = " + time + " x = " + x[0]);
			time = rkg(time, h, x, dx, g, n, snx);
		}
	}
			// 微係数の計算
	int snx(double time, double[] x, double[] dx)
	{
		dx[0] = x[1];
		dx[1] = -2.0 * x[0] - 3.0 * x[1] + 1.0;
		return 0;
	}

	/********************************************/
	/* ルンゲ・クッタ法  dx/dt=f(t,x)           */
	/*      time : 現在の時間                   */
	/*      h : 時間刻み幅                      */
	/*      x : 現在の状態                      */
	/*      dx : 微係数（f(t,x):snxで計算する） */
	/*      g : 作業域（g[4][n]）               */
	/*      n : 微分方程式の次数                */
	/*      fn : 関数値を計算する関数           */
	/*      return : time+h                     */
	/********************************************/
	double rkg(double time, double h, double[] x, double[] dx, double[][] g,
	           int n, Func<double, double[], double[], int> fn)
	{
		double h2 = 0.5 * h;

		fn(time, x, dx);
		for (int i1 = 0; i1 < n; i1++)
			g[0][i1] = h * dx[i1];

		time += h2;
		for (int i1 = 0; i1 < n; i1++)
			g[1][i1] = x[i1] + 0.5 * g[0][i1];
		fn(time, g[1], dx);
		for (int i1 = 0; i1 < n; i1++)
			g[1][i1] = h * dx[i1];

		for (int i1 = 0; i1 < n; i1++)
			g[2][i1] = x[i1] + 0.5 * g[1][i1];
		fn(time, g[2], dx);
		for (int i1 = 0; i1 < n; i1++)
			g[2][i1] = h * dx[i1];

		time += h2;
		for (int i1 = 0; i1 < n; i1++)
			g[3][i1] = x[i1] + g[2][i1];
		fn(time, g[3], dx);
		for (int i1 = 0; i1 < n; i1++)
			g[3][i1] = h * dx[i1];

		for (int i1 = 0; i1 < n; i1++)
			x[i1] = x[i1] + (g[0][i1] + 2.0 * g[1][i1] + 2.0 * g[2][i1] + g[3][i1]) / 6.0;

		return time;
	}
}

''''''''''''''''''''''''''''''''''
' 微分方程式（ルンゲ・クッタ法） '
'      例：d2x/dt+3dx/dt+2x=1    '
'      coded by Y.Suganuma       '
''''''''''''''''''''''''''''''''''
Module Test

	Sub Main()
				'
				' 初期設定
				'
		Dim n As Integer   = 2
		Dim time As Double = 0.0
		Dim h As Double    = 0.01
		Dim x() As Double  = {0.0, 0.0}
		Dim dx(n) As Double
		Dim g0(2) As Double
		Dim g1(2) As Double
		Dim g2(2) As Double
		Dim g3(2) As Double
				' 微係数の計算
		Dim snx = Function(tm, v, dv)
			dv(0) = v(1)
			dv(1) = -2.0 * v(0) - 3.0 * v(1) + 1.0
			Return 0
		End Function
				'
				' 計算と出力
				'
		For i1 As Integer = 0 To 100
			Console.WriteLine("time = " & time & " x = " & x(0))
			time = rkg(time, h, x, dx, g0, g1, g2, g3, n, snx)
		Next

	End Sub

	''''''''''''''''''''''''''''''''''''''''''''
	' ルンゲ・クッタ法  dx/dt=f(t,x)           '
	'      time : 現在の時間                   '
	'      h : 時間刻み幅                      '
	'      x : 現在の状態                      '
	'      dx : 微係数（f(t,x):snxで計算する） '
	'      g0,g1,g2,g3 : 作業域（gi(n)）       '
	'      n : 微分方程式の次数                '
	'      fn : 関数値を計算する関数           '
	'      return : time+h                     '
	''''''''''''''''''''''''''''''''''''''''''''
	Function rkg(time As Double, h As Double, x() As Double, dx() As Double,
	             g0() As Double, g1() As Double, g2() As Double, g3() As Double,
	             n As Integer, fn As Func(Of Double, Double(), Double(), Integer))

		Dim h2 As Double = 0.5 * h

		fn(time, x, dx)
		For i1 As Integer = 0 To n-1
			g0(i1) = h * dx(i1)
		Next

		time += h2
		For i1 As Integer = 0 To n-1
			g1(i1) = x(i1) + 0.5 * g0(i1)
		Next
		fn(time, g1, dx)
		For i1 As Integer = 0 To n-1
			g1(i1) = h * dx(i1)
		Next

		For i1 As Integer = 0 To n-1
			g2(i1) = x(i1) + 0.5 * g1(i1)
		Next
		fn(time, g2, dx)
		For i1 As Integer = 0 To n-1
			g2(i1) = h * dx(i1)
		Next

		time += h2
		For i1 As Integer = 0 To n-1
			g3(i1) = x(i1) + g2(i1)
		Next
		fn(time, g3, dx)
		For i1 As Integer = 0 To n-1
			g3(i1) = h * dx(i1)
		Next

		For i1 As Integer = 0 To n-1
			x(i1) = x(i1) + (g0(i1) + 2.0 * g1(i1) + 2.0 * g2(i1) + g3(i1)) / 6.0
		Next

		Return time

	End Function

End Module