global m n C L R tt b

m = 7; n = 3;
C = 0.2; L = 0.15; R = 0.005;
b = [0;1;-i]/2/L;
N = 512;

tau = input('tau = ');
tt  = [0, 2*pi/(m-1)*[1/2:(m-3/2)]+tau, 2*pi];

f0 = solution(zeros(n,1));
F = [];
for x = eye(n),
   F = [F, f0 - solution(x)];
end;
x0 = F\f0;
[ff, xx, uu, M] = solution(x0);

delta = 2*pi/N; 
k = 0; xtab = [];
for j = 0:N-1, 
   t = j*delta;
   if (tt(k+1) <= t),
      k = k+1;
      A = matrix(k-1);
      aux = expm(A*(t-tt(k))) ...
            *(xx(:,k)+2*real(uu(:,k)*exp(i*tt(k))));
      E = expm(A*delta);
   else
      aux = E*aux;
   end;
   xtab = [xtab, aux-2*real(uu(:,k)*exp(i*t))];
end;

figure(1); plot(xtab(1,:)');
figure(2); plot(xtab(2:3,:)');

cc = fft(1/N*xtab'); cc(1:32,:)

