%  Solving the linear diferential equation  X' = -i A(t) X
% using the 8th-order MAGNUS integrators from Blanes et al, BIT (2000)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global n
% Dimension of the matrix
n=20;
% Time of integration: t0 (initial time), tf (final time)
t0 = -2; tf=t0+4*2*pi;
% Evaluations
ipart = round(tf-t0)*300;
for ii = 1:1;
% Initial conditions
  xm8 = eye(n);
% Time-step
  h = (tf-t0)/ipart;
% Quadrature points
  v1 = sqrt((3+2*sqrt(6/5))/7)/2;   v2 = sqrt((3-2*sqrt(6/5))/7)/2;
  c1 = 1/2 - v1;   c2 = 1/2 - v2;   c3 = 1/2 + v2;   c4 = 1/2 + v1;
  t = t0;
  for j2 = 1:ipart;
%  Evaluation of the matrices at the Gaussian quadratures
  t1 = t + c1*h;   t2 = t + c2*h;   t3 = t + c3*h;   t4 = t + c4*h;
  A1 = -i*H(t1);   A2 = -i*H(t2);   A3 = -i*H(t3);   A4 = -i*H(t4);
    
   ww1 = 1/2 - sqrt(5/6)/6;
   ww2 = 1/2 + sqrt(5/6)/6;
   S1 = A1 + A4;
   S2 = A2 + A3;
   R1 = A4 - A1;
   R2 = A3 - A2;
  B0 = h*(   ww1*S1   +    ww2*S2)/2;
  B1 = h*(v1*ww1*R1   + v2*ww2*R2)/2;
  B2 = h*(v1^2*ww1*S1 + v2^2*ww2*S2)/2;
  B3 = h*(v1^3*ww1*R1 + v2^3*ww2*R2)/2;
   SU1 = -B0*38/5 + 24*B2;
  Q1 = SU1*B3 - B3*SU1;
   SU2 = B0*63/5 - 84*B2;
   SU3 = -B1*5/28 + B3;
  Q2 = SU2*SU3 - SU3*SU2;
   SU4 = B0*19/28 - B2*15/7;
   SU5 = B2 + Q1*61/588 - Q2/12;
   CO1 = B0*SU5 - SU5*B0;
  Q3 = SU4*CO1 - CO1*SU4;
   SU6 = Q1*20/7 + 10*Q2;
  Q4 = B3*SU6 - SU6*B3;
   CO2 = B2*Q1 - Q1*B2;
   SU7 = -B0*6025/4116 + B2*2875/343;
  Q5 = SU7*CO2 - CO2*SU7;
   SU8 = (Q3 + Q4)*20/7 + Q5*820/189;
  Q6 = B3*SU8 - SU8*B3;
   SU9 = -(Q3 - Q4/3 + Q5)/42;
   CO3 = B0*SU9 - SU9*B0;
  Q7 = B0*CO3 - CO3*B0;
  M8 = B0 + Q1 + Q2 + Q3 + Q4 + Q5 + Q6 + Q7 ;
% Exponential
  y0 =expm(M8)*xm8;
    xm8 = y0;
    t = t +h;
  end
xex=xm8
end
