--- a/+time/+rk/rungekutta.m	Mon Nov 26 16:08:54 2018 -0800
+++ b/+time/+rk/rungekutta.m	Mon Nov 26 16:12:27 2018 -0800
@@ -2,15 +2,18 @@
 % starting from @arg v and where the function F(v,t) gives the
 % time derivatives. coeffs is a struct holding the RK coefficients
 % for the specific method.
-function v = rungekutta(v, t , dt, F, coeffs)
+% Also returns the stage approximations (V) and stage rates (K).
+function [v, V, K] = rungekutta(v, t , dt, F, coeffs)
     % Compute the intermediate stages k
-    k = zeros(length(v), coeffs.s);
+    K = zeros(length(v), coeffs.s);
+    V = zeros(length(v), coeffs.s);
     for i = 1:coeffs.s
         u = v;
         for j = 1:i-1
-            u = u + dt*coeffs.a(i,j)*k(:,j);
+            u = u + dt*coeffs.a(i,j)*K(:,j);
         end
-        k(:,i) = F(u,t+coeffs.c(i)*dt);
+        V(:,i) = u;
+        K(:,i) = F(u,t+coeffs.c(i)*dt);
     end
     % Compute the updated solution as a linear combination
     % of the intermediate stages.