function dydt = sir_age_model(t, y, parameters)
    % Unpack parameters
    b = parameters(1);
    c_11 = parameters(2);
    c_12 = parameters(3);
    c_21 = parameters(4);
    c_22 = parameters(5);
    gamma = parameters(6);

    % Unpack state variables
    S1 = y(1);
    I1 = y(2);
    R1 = y(3);
    S2 = y(4);
    I2 = y(5);
    R2 = y(6);

    % Total population in each group
    N1 = S1 + I1 + R1;
    N2 = S2 + I2 + R2;

    % Force of infection
    lambda_1 = b * (c_11 * I1 / N1 + c_12 * I2 / N2);
    lambda_2 = b * (c_21 * I1 / N1 + c_22 * I2 / N2);

    % Differential equations
    dS1 = -lambda_1 * S1;
    dI1 = lambda_1 * S1 - gamma * I1;
    dR1 = gamma * I1;
    dS2 = -lambda_2 * S2;
    dI2 = lambda_2 * S2 - gamma * I2;
    dR2 = gamma * I2;

    % Output derivative
    dydt = [dS1; dI1; dR1; dS2; dI2; dR2];
end