function [ dSdt ] = Lsquareddynamics(t,S)
% this function 
global r1 rn1 r2 rn2 r3 

dSdt = zeros(size(S));

a = S(1); %
b = S(2); % 
c = S(3); % 

dadt  = -r1*a + rn1*b - r3*a*c;
dbdt  = r1*a - rn1*b + rn2*c - r2*b + 2*r3*a*c;
dcdt  = -rn2*c + r2*b - r3*a*c;

dSdt(1) = dadt;
dSdt(2) = dbdt;
dSdt(3) = dcdt;
end