function dX = systemcl(t, X)
% Equation describing the evolution of the CL model.
% The model evolves according to:   
%     x'' = Ax
% Where x is the mass position vector and A is the
% matrix describing the spring system.
% Using (u' = x'') to transform the 2nd order equation into 
% two first order equations, we obtain:
%     dx = u, u' = Ax 
%     dX  = [x' u']
persistent A M P N;

if isempty(A)
    M = size(X,1);
    P = (M/2); N = P - 1;
    A = gimmeA(M,P,N);    
end

dX = zeros(M,1); 
dX(1:P)= X(P+1:end);     % x' = u
dX(P+1:end) = A*X(1:P);  % u' = Ax

function A = gimmeA(M,P,N)
    omega = 1; gamma = 1; m = 1;  w = 10/N; g = sqrt(40*gamma/(N*pi));
    K = N * (g^2) / m;
    A = eye(P).*(repmat(-w*[0:N].^2,P,1));
    A(1,:) = [-((omega^2)+K/m)  (g/m)*w*[1:N]];
    A(2:end,1) = ((g/m)*w*[1:N])';