function res = muschel_f(t, y)

global N s D_a D_b D_c b_a b_b b_c s_a s_b s_c r_a r_b r_c Lap indA indB indC

a = y(indA);
b = y(indB);
c = y(indC);

aq = s .* b .* ((a.^2 + b_a) ./ (1 + s_a * a.^2)) ./ (s_b + s_c * c);

olddecaydiffA = -r_a * a + D_a * (Lap * a);
olddecaydiffB = -r_b * b + D_b * (Lap * b);
olddecaydiffC = -r_c * c + D_c * (Lap * c);

res = [ olddecaydiffA + aq; olddecaydiffB - aq + b_b; olddecaydiffC + r_c * a ];