%  Solving the linear diferential equation  X' = i A(t) X
% taking  X = q + i p
% using different  split-MAGNUS methods of ORDER 4 and 6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global n
% n: Dimension of the matrix (defined in Ex.m)
% Time of integration: t0 (initial time), tf (final time)
t0 = -2; tf=t0+4*2*pi;
% Number of periods
per = tf-t0;
% Quadrature points: 6th-order GL
c = [1/2 - sqrt(3/20), 1/2, 1/2 + sqrt(3/20)];
w = [5,8,5]/18;
c1 = c(1);     c2 = c(2);    c3 = c(3);

imet=[1,2,3];
for mm=1:3;
y1 = [];  ww = [];
 [a,b, nex]=coefsSplit(imet(mm));
% Evaluations per period
  ipart = 2;
%%%%%%%  repeated integrations for different time steps  %%%%%%%%%%%%%
 for ip=1:8
% Initial conditions
  xm6 = eye(n);
  q=real(xm6);       p=imag(xm6);
  t = t0;
  % Time-step
  h = (tf-t0)/ipart;
   for j2 = 1:ipart;

%  Evaluation of the matrices at the 6th-order Gaussian quadratures
    t1 = t + c1*h;    t2 = t + c2*h;    t3 = t + c3*h;
%%%%%%%%%%%%%%%%%%%%%
    H1 = H(t1);    H2 = H(t2);    H3 = H(t3);
    for i2=1:nex+1
      q=q+h*(a(i2,1)*H1+a(i2,2)*H2+a(i2,3)*H3)*p;  
      p=p-h*(b(i2,1)*H1+b(i2,2)*H2+b(i2,3)*H3)*q; 
    end
    
    xm6=q + p*1i;
    t = t +h;
   end
   
   error = log10(norm(xm6 - xex));
   y1 = [y1 error];
   ww = [ww log10(ipart*nex)];
   ipart=ipart*2;
 end
 
% Errors versus cost
  [ww',y1']
% Plots of the efficiency figures
  hold on
  if (mm==1)
   plot(ww,y1,'k-o','LineWidth',2.)
  elseif (mm==2)
   plot(ww,y1,'k-s','LineWidth',2.)
  else 
   plot(ww,y1,'k-p','LineWidth',2.)
  end
end
axis([1.9 3.5 -8 1])
box
xlabel('{\fontsize{18} log_{10}(N. products)}')
ylabel('{\fontsize{16}log_{10}(error)}')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %  figure(1)
% title('{\fontsize{14}\epsilon=0,  T_0=10}')
% % %  figure(2)
% title('{\fontsize{14}\epsilon=1/10,  \delta=1/10,  T_0=5}')
% % %  figure(3)
% title('{\fontsize{14}\epsilon=1/2,  \delta=1,  T_0=5}')
% % %  figure(4)
% title('{\fontsize{14}\epsilon=2,  \delta=5,  T_0=1}')

