% The matrix of flow equations for the gaussian flow 
% the components are y(1)=Kc(l) y(2) = vc(l) y(3)=y(l) y(4)=eta(l) y(5)=lambda(l) 
% y(6) = D(k1) y(7) = D(k2) ...
function dA = gausFullMatrix(l,A1)
    A = real(A1);
    Kc = A(1); vc = A(2); Ks = A(3); vs = A(4); y = A(5); eta = A(6); lam = A(7);
    Db = A(8:length(A)); 
    N = length(A)-7;

    k = (2*pi/N).*(0:N-1)';
    r1 = (vs/vc)^(2-Kc);
    r2 = (vs/vc)^(-Kc);
    r3 = 1+(1/Kc^2);
    r4 = 1-(1/Kc^2);   
    r3s = 1+(1/Ks^2);
    r4s= 1-(1/Ks^2);
    f  = -0.25*(r1*r4*Kc+r2*r4s*Ks);
    
    dKc = -0.5*(vc/vs)*(Kc^2)*sum(Db.*(1-0.5*r1*r3*(3-Kc-Ks-y-k.^2)./(1-k.^2)));
    dvc = -0.5*(Kc/vs)*(vc^2)*sum(Db.*(1+0.5*r1*r4*(3-Kc-Ks-y-k.^2)./(1-k.^2)));
    dKs = -0.5*Ks^2*(y^2+sum(Db.*(1-0.5*r2*(r3s.*(3-Kc-Ks-y-k.^2)+4*y)./(1-k.^2))));
    dvs = -0.5*Ks*vs*sum(Db.*(1+0.5*r2*r4s*(3-Kc-Ks-y-k.^2)./(1-k.^2)));
    dy  = (2-2*Ks)*y - sum(Db.*(1-r2*(1+0.5*r2*r4s*(3-Kc-Ks-y-k.^2)./(1-k.^2)-(r3s*y*Ks^2-2*(Kc-Ks+y))./(2-2*k.^2))));
    dDb = (3-Kc-Ks-y-k.^2).*Db;
    dlam= -lam*eta*sum(Db.*k.^2.*(1-f-f*(2-Kc-Ks-y-k.^2).*(1-exp(-l*(1-k.^2)))./(1-k.^2)));
    deta = eta*dlam/lam;
    
    par = [l Kc Db(1)]
    %if r5==0
    %    dlmwrite('r5Log.m',[l min(abs(r5)) max(abs(r5)) lambda Db(1) Kc],'-append','roffset',1,'delimiter','\t');
    %end
    
dA = vertcat([dKc; dvc; dKs; dvs; dy; deta; dlam],dDb);
