--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/+time/rungekuttaRV.m	Thu Jan 17 10:44:12 2019 +0100
@@ -0,0 +1,30 @@
+% Takes one time step of size dt using the rungekutta method
+% starting from v and where the function F(v,t,RV) gives the
+% time derivatives. coeffs is a struct holding the RK coefficients
+% for the specific method. RV is the residual viscosity which is updated
+% in between the stages and after the updated solution is computed.
+function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs)
+    % Move one stage outside to avoid branching for updating the
+    % residual inside the loop.
+    k = zeros(length(v), coeffs.s);
+    k(:,1) = F(v,t,RV.evaluate(v,DvDt(v)));
+
+    % Compute the intermediate stages k
+    for i = 2:coeffs.s
+        u = v;
+        for j = 1:i-1
+            u = u + dt*coeffs.a(i,j)*k(:,j);
+        end
+        %RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization
+        k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluate(u,DvDt(u)));
+    end
+
+    % Compute the updated solution as a linear combination
+    % of the intermediate stages.
+    u = v;
+    for i = 1:coeffs.s
+        u = u + dt*coeffs.b(i)*k(:,i);
+    end
+    %RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization
+    v = u;
+end