% 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 , matrixexp ,tol)



if isa(D,'function_handle')
    c1 = 1/2 - sqrt(3)/6;
    c2 = 1/2 + sqrt(3)/6;
    
    A1 = D(t +c1*k);
    A2 = D(t + c2*k);
    Omega = 1/2*(A1 + A2) + sqrt(3)*k/12*(A1*A2-A2*A1);
else
    Omega = D;
end


switch matrixexp
    case 'expm'
        v = expm(k*Omega)*v;
    case 'Arnoldi'
        v = time.expint.expm_Arnoldi(-Omega,v,k,tol,100);
    otherwise
        error('No such matrix exponential evaluation')
        
end
end