function res = muschel_f_sp76(t, y)

global N s D_a D_b D_c D_d D_e D_f
global b_a b_b b_c b_d b_e b_f
global c_a c_b c_c c_d c_e c_f
global r_a r_b r_c r_d r_e r_f
global s_a s_b s_c s_d s_e s_f
global Lap indA indB indC indD indE indF

a = y(indA);
b = y(indB);
c = y(indC);
d = y(indD);
e = y(indE);
f = y(indF);

a_star = a.^2 ./ (1 + s_a * a.^2) + b_a;
e_star = e.^2 ./ (1 + s_e * e.^2) + b_e;

na = s .* b .* a_star -r_a * a -c_e .* e .* a + (D_a*Lap*a);
nb = b_b * (1 + s_b * c + c_b * d) - s .* b .* a_star -r_b * b + (D_b*Lap*b);
nc = r_c * a .* ((c.^2 + b_c) ./ (s_d + d)) - r_c * c + (D_c*Lap*c);
nd = r_c * a .* (c.^2 + b_c) - r_d * d + b_d + (D_d*Lap*d);
ne = r_e * f .* e_star - r_e * e + s_f * a + (D_e*Lap*e);
nf = b_f * a + c_f - r_e * f .* e_star - r_f * f + (D_f*Lap*f);

res = [ na; nb; nc; nd; ne; nf ];