% 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