function dY = brusselator(time,Y,FLAG,struct)

%space allocation for the derivative
%(first N derivatives for [X], last N derivatives for [Y])
dY=zeros(struct.N.*2,1);

%only one cell?
if struct.N==1
    %yes
    %Y(1) is [X]
    dY(1)=struct.k1A+struct.k2*Y(1)^2*Y(2)-struct.k3B*Y(1)-struct.k4*Y(1);
    %Y(2) is [Y]
    dY(2)=-struct.k2*Y(1)^2*Y(2)+struct.k3B*Y(1);
else
    %dynamic equations for [X] in the N cells
    for i=1:1:struct.N
        %is i an "edge" cells?
        if i==1
            %i is the initial cell in the line (no left neighbour)
            dY(i)=struct.k1A+struct.k2*Y(i)^2*Y(i+struct.N)-struct.k3B*Y(i)-struct.k4*Y(i)+struct.Dx*(-2*Y(i)+Y(i+1));
        elseif i==struct.N
            %i is the final cell in the line (no right neighbour)
            dY(i)=struct.k1A+struct.k2*Y(i)^2*Y(i+struct.N)-struct.k3B*Y(i)-struct.k4*Y(i)+struct.Dx*(Y(i-1)-2*Y(i));
        else
            dY(i)=struct.k1A+struct.k2*Y(i)^2*Y(i+struct.N)-struct.k3B*Y(i)-struct.k4*Y(i)+struct.Dx*(Y(i-1)-2*Y(i)+Y(i+1));
        end
    end
    %dynamic equations for [Y] in the N cells
    for i=1:1:struct.N
        %is i an "edge" cells?
        if i==1
            %i is the initial cell in the line (no left neighbour)
            dY(i+struct.N)=-struct.k2*Y(i)^2*Y(i+struct.N)+struct.k3B*Y(i)+struct.Dy*(-2*Y(i+struct.N)+Y(i+1+struct.N));
        elseif i==struct.N
            %i is the final cell in the line (no right neighbour)
            dY(i+struct.N)=-struct.k2*Y(i)^2*Y(i+struct.N)+struct.k3B*Y(i)+struct.Dy*(Y(i-1+struct.N)-2*Y(i+struct.N));
        else
            dY(i+struct.N)=-struct.k2*Y(i)^2*Y(i+struct.N)+struct.k3B*Y(i)+struct.Dy*(Y(i-1+struct.N)-2*Y(i+struct.N)+Y(i+1+struct.N));
        end
    end
end