From aa3a086cdab9860619e3b271a29fa5a1a803ff10 Mon Sep 17 00:00:00 2001
From: jgoppert <james.goppert@gmail.com>
Date: Thu, 12 Nov 2015 15:13:17 -0500
Subject: [PATCH] Work on rk4 interface.

---
 matrix/integration.hpp | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/matrix/integration.hpp b/matrix/integration.hpp
index 1979bd8753..67eae7a453 100644
--- a/matrix/integration.hpp
+++ b/matrix/integration.hpp
@@ -4,10 +4,11 @@
 
 namespace matrix {
 
-template<typename Type, size_t M>
+template<typename Type, size_t M, size_t N>
 int integrate_rk4(
-    Vector<Type, M> (*f)(Type, Vector<Type, M>),
+    Vector<Type, M> (*f)(Type, const Vector<Type, M> &x, const Vector<Type, N> & u),
     const Vector<Type, M> & y0,
+    const Vector<Type, N> & u,
     Type t0,
     Type h,
     Vector<Type, M> & y1
@@ -15,10 +16,10 @@ int integrate_rk4(
 {
     // https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
     Vector<Type, M> k1, k2, k3, k4;
-    k1 = f(t0, y0);
-    k2 = f(t0 + h/2, y0 + k1*h/2);
-    k3 = f(t0 + h/2, y0 + k2*h/2);
-    k4 = f(t0 + h, y0 + k3*h);
+    k1 = f(t0, y0, u);
+    k2 = f(t0 + h/2, y0 + k1*h/2, u);
+    k3 = f(t0 + h/2, y0 + k2*h/2, u);
+    k4 = f(t0 + h, y0 + k3*h, u);
     y1 = y0 + (k1 + k2*2 + k3*2 + k4)*(h/6);
     return 0;
 }