noisy_deriv.m

function dydt = noisy_deriv(y, nu, mu, epsilon, beta, gamma, Gamma, q)
    S = y(1);
    I = y(2);
    R = y(3);
    C = y(4);
    noise = 25 * randn(4);
    dSdt = nu - (beta*I + epsilon*beta*C)*S - mu*S + noise(1);
    dIdt = (beta*I + epsilon*beta*C)*S - gamma*I - mu*I + noise(2);
    dRdt = gamma*(1-q)*I + Gamma*C - mu*R + noise(3);
    dCdt = gamma*q*I - Gamma*C - mu*C + noise(4);
    dydt = [dSdt; dIdt; dRdt; dCdt];
end