diff -r 878652b22157 -r 679f4ddd982f +time/Rungekutta.m
--- a/+time/Rungekutta.m	Mon Nov 26 16:12:27 2018 -0800
+++ b/+time/Rungekutta.m	Mon Nov 26 16:23:27 2018 -0800
@@ -1,45 +1,58 @@
 classdef Rungekutta < time.Timestepper
     properties
         F       % RHS of the ODE
-        k       % Time step
+        dt      % Time step
         t       % Time point
         v       % Solution vector
         n       % Time level
         scheme  % The scheme used for the time stepping, e.g rk4, rk6 etc.
+        coeffs  % Butcher tableau coefficients
+        V       % All stage approximations in most recent time step
+        K       % All stage rates in most recent time step
     end
 
 
     methods
         % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with
-        % timestep k and initial conditions v = v0
-        function obj = Rungekutta(F, k, t0, v0, method)
+        % timestep dt and initial conditions v = v0
+        function obj = Rungekutta(F, dt, t0, v0, method)
             default_arg('method',"rk4");
             obj.F = F;
-            obj.k = k;
+            obj.dt = dt;
             obj.t = t0;
             obj.v = v0;
             obj.n = 0;
+
+            % Extract the coefficients for the specified method
+            % used for the RK updates from the Butcher tableua.
+            [s,a,b,c] = time.rk.butcherTableau(method);
+            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+
             % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation
             % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
             if (method == "rk4")
                 obj.scheme = @time.rk.rungekutta_4;
             else
-                % Extract the coefficients for the specified method
-                % used for the RK updates from the Butcher tableua.
-                [s,a,b,c] = time.rk.butcherTableau(method);
-                coeffs = struct('s',s,'a',a,'b',b,'c',c);
-                obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
+                obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, obj.coeffs);
             end
         end
 
-        function [v,t] = getV(obj)
+        % v: Current solution
+        % t: Current time
+        % V: All stage approximations in most recent time step
+        % K: All stage rates in most recent time step
+        % T: Time points (corresponding to V and K) in most recent time step
+        function [v,t,V,T,K] = getV(obj)
             v = obj.v;
             t = obj.t;
+            V = obj.V;
+            K = obj.K;
+            T = obj.t + obj.dt*obj.coeffs.b;
         end
 
         function obj = step(obj)
-            obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
+            [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.dt, obj.F);
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end