--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+expint/Magnus_4.m	Fri Jun 02 17:46:01 2017 +0200
@@ -0,0 +1,20 @@
+% Takes one time step of size k using a fourth order magnus integrator
+% starting from v_0 and where the function F(v,t) gives the
+% time derivatives.
+function v = Magnus_4(v,D, t , k)
+
+if isa(D,'function_handle')
+   % v = krylov(k*D(t +k/2*t),v);
+   c1 = 1/2 - sqrt(3)/6;
+   c2 = 1/2 + sqrt(3)/6;  
+   
+   A1 = D(t +c1*k);
+   A2 = D(t + c2*k);
+   Omega = k/2*(A1 + A2) + sqrt(3)*k^2/12*(A1*A2-A2*A1);
+   v = expm(Omega)*v;
+else
+   %v = krylov(k*D,v);
+   v = expm(k*D)*v;
+end
+
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/Magnus4.m	Fri Jun 02 17:46:01 2017 +0200
@@ -0,0 +1,43 @@
+classdef Magnus4 < time.Timestepper
+    properties
+        D
+        S
+        F
+        k
+        t
+        v
+        m
+        n
+    end
+
+
+    methods
+        function obj = Magnus4(D, k, t0, v0)
+            obj.D = D;
+            obj.k = k;
+            obj.t = t0;
+            obj.v = v0;
+            obj.m = length(v0);
+            obj.n = 0;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            obj.v = time.expint.Magnus_4(obj.v,obj.D, obj.t, obj.k);
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+        end
+    end
+
+
+    methods (Static)
+        function k = getTimeStep(lambda)
+            k = rk4.get_rk4_time_step(lambda);
+        end
+    end
+
+end
\ No newline at end of file