diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index b6004e90..6369d13e 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -112,6 +112,7 @@
+
diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs
new file mode 100644
index 00000000..dd3ab6fe
--- /dev/null
+++ b/src/Numerics/OdeSolvers/OdeSolvers.cs
@@ -0,0 +1,156 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using MathNet.Numerics.LinearAlgebra;
+using System;
+
+namespace MathNet.Numerics.OdeSolvers
+{
+ ///
+ /// ODE Solver Algorithms
+ ///
+ public static class RungeKutta
+ {
+ ///
+ /// Second Order Runge-Kutta method
+ ///
+ /// initial value
+ /// start time
+ /// end time
+ /// Number of subintervals
+ /// ode function
+ /// approximations
+ public static double[] SecondOrder(double y0, double start, double end, int N, Func f)
+ {
+ double dt = (end - start) / (N - 1);
+ double k1 = 0;
+ double k2 = 0;
+ double t = start;
+ double[] y = new double[N];
+ y[0] = y0;
+ for (int i = 1; i < N; i++)
+ {
+ k1 = f(t, y0);
+ k2 = f(t + dt, y0 + k1 * dt);
+ y[i] = y0 + dt * 0.5 * (k1 + k2);
+ t += dt;
+ y0 = y[i];
+ }
+ return y;
+ }
+
+ ///
+ /// Fourth Order Runge-Kutta method
+ ///
+ /// initial value
+ /// start time
+ /// end time
+ /// number of subintervals
+ /// ode function
+ /// approximations
+ public static double[] FourthOrder(double y0, double start, double end, int N, Func f)
+ {
+ double dt = (end - start) / (N - 1);
+ double k1 = 0;
+ double k2 = 0;
+ double k3 = 0;
+ double k4 = 0;
+ double t = start;
+ double[] y = new double[N];
+ y[0] = y0;
+ for (int i = 1; i < N; i++)
+ {
+ k1 = f(t, y0);
+ k2 = f(t + dt / 2, y0 + k1 * dt / 2);
+ k3 = f(t + dt / 2, y0 + k2 * dt / 2);
+ k4 = f(t + dt, y0 + k3 * dt);
+ y[i] = y0 + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
+ t += dt;
+ y0 = y[i];
+ }
+ return y;
+ }
+
+ ///
+ /// Second Order Runge-Kutta to solve ODE SYSTEM
+ ///
+ /// initial vector
+ /// start time
+ /// end time
+ /// number of subintervals
+ /// ode function
+ /// approximations
+ public static Vector[] SecondOrder(Vector y0, double start, double end, int N, Func, Vector> f)
+ {
+ double dt = (end - start) / (N - 1);
+ Vector k1, k2;
+ Vector[] y = new Vector[N];
+ double t = start;
+ y[0] = y0;
+ for (int i = 1; i < N; i++)
+ {
+ k1 = f(t, y0);
+ k2 = f(t, y0 + k1 + dt);
+ y[i] = y0 + 0.5 * (k1 + k2);
+ t += dt;
+ y0 = y[i];
+ }
+ return y;
+ }
+
+ ///
+ /// Fourth Order Runge-Kutta to solve ODE SYSTEM
+ ///
+ /// initial vector
+ /// start time
+ /// end time
+ /// number of subintervals
+ /// ode function
+ /// approximations
+ public static Vector[] FourthOrder(Vector y0, double start, double end, int N, Func, Vector> f)
+ {
+ double dt = (end - start) / (N - 1);
+ Vector k1, k2, k3, k4;
+ Vector[] y = new Vector[N];
+ double t = start;
+ y[0] = y0;
+ for (int i = 1; i < N; i++)
+ {
+ k1 = f(t, y0);
+ k2 = f(t + dt / 2, y0 + k1 * dt / 2);
+ k3 = f(t + dt / 2, y0 + k2 * dt / 2);
+ k4 = f(t + dt, y0 + k3 * dt);
+ y[i] = y0 + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
+ t += dt;
+ y0 = y[i];
+ }
+ return y;
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs
new file mode 100644
index 00000000..d8b3542b
--- /dev/null
+++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs
@@ -0,0 +1,91 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using NUnit.Framework;
+using System;
+using MathNet.Numerics.OdeSolvers;
+using System.Linq;
+
+namespace MathNet.Numerics.UnitTests.OdeSolvers
+{
+ ///
+ /// ODE Solver tests.
+ ///
+ [TestFixture, Category("OdeSolver")]
+ public class OdeSolverTest
+ {
+ ///
+ /// Runge-Kutta second order method for first order ODE.
+ ///
+ [Test]
+ public void RK2Test()
+ {
+ Func ode = (t, y) => t + 2 * y * t;
+ Func sol = (t) => 0.5 * (Math.Exp(t * t) - 1);
+ double ratio = double.NaN;
+ double error = 0;
+ double oldError = 0;
+ for (int k = 0; k < 4; k++)
+ {
+ double y0 = 0;
+ double[] y_t = RungeKutta.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode);
+ error = Math.Abs(sol(2) - y_t.Last());
+ if (oldError != 0)
+ ratio = Math.Log(oldError / error, 2);
+ oldError = error;
+ Console.WriteLine(string.Format("{0}, {1}", error, ratio));
+ }
+ Assert.AreEqual(2, ratio, 0.01);// Check error convergence order
+ }
+
+ ///
+ /// Runge-Kutta fourth order method for first order ODE.
+ ///
+ [Test]
+ public void RK4Test()
+ {
+ Func ode = (t, y) => t + 2 * y * t;
+ Func sol = (t) => 0.5 * (Math.Exp(t * t) - 1);
+ double ratio = double.NaN;
+ double error = 0;
+ double oldError = 0;
+ for (int k = 0; k < 4; k++)
+ {
+ double y0 = 0;
+ double[] y_t = RungeKutta.FourthOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode);
+ error = Math.Abs(sol(2) - y_t.Last());
+ if (oldError != 0)
+ ratio = Math.Log(oldError / error, 2);
+ oldError = error;
+ Console.WriteLine(string.Format("{0}, {1}", error, ratio));
+ }
+ Assert.AreEqual(4, ratio, 0.01);// Check error convergence order
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 735d4a80..e07050ea 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -362,6 +362,7 @@
+