See

<project-file type=“source”/> <content> #include “equations.h” #include <gsl/gsl_errno.h> #include “ke_index.h”

This function defines the index of the variable when its degree is k, and the infectious degree (we used to call it the effective degree) is e. /* inline int ke_index(int k,int e) { return 1)

1)
k-1)*(k-1+3)/2 + e + 1); } */ inline double V(const double X[], int n, int i, int max) {
  if (n < 0) return 0;
  if (n > max) return 0;
  if (i > n) return 0;
  if (i < 0) return 0;
  return X[ke_index(n, i)];
} int ode(double t, const double y[], double dydt[], void *params) {
  const ODEParameters *p = (ODEParameters*) params;
  int Mh = p->MaxHouseholdSize;
  int mh = Mh * (Mh + 3) / 2 + 1;
  int Msi = p->MaxExternalDegree;
  int msi = Msi * (Msi + 3) / 2 + 1;
  int m = mh + msi * 2;
  const double *H = y, *S = &y[mh], *I = &y[mh + msi];
  double beta_h = p->BetaHousehold;
  double beta_e = p->BetaExternal;
  double gamma = p->Gamma;
  double external_infection_rate = 0;
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++)
          external_infection_rate += beta_e * i * V(S, k, i, Msi);
  double household_susceptible = 0;
  for (int n = 0; n <= Mh; n++)
      for (int j = 0; j <= n; j++)
          household_susceptible += pnj(n, j) * V(H, n, j, Mh);
  // the H equations
  for (int n = 0; n <= Mh; n++)
      for (int j = 0; j <= n; j++) {
          dydt[ke_index(n, j)] =
                  beta_h * ((n - j + 1) * (j - 1) * V(H, n, j - 1, Mh) - (n - j) * j * V(H, n, j, Mh)) +
                  gamma * ((j + 1) * V(H, n + 1, j + 1, Mh) - j * V(H, n, j, Mh)) +
                  external_infection_rate / household_susceptible *
                      (pnj(n, j - 1) * V(H, n, j - 1, Mh)  - pnj(n, j) * V(H, n, j, Mh));
      }
  double incidence_s_neighbors = 0;
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          incidence_s_neighbors += beta_e * (k - i) * i * V(S, k, i, Msi);
      }
  double susceptible_s_neighbors = 0;
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          susceptible_s_neighbors += (k - i) * V(S, k, i, Msi);
      }
  double internal_infection_rate = 0;
  for (int n = 0; n <= Mh; n++)
      for (int j = 0; j <= n; j++)
          internal_infection_rate += beta_h * (n - j) * j * V(H, n, j, Mh);
  double all_s = 0;
  for (int j = 0; j < msi; j++) all_s += S[j];
  // The S equations
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          int s = k - i;
          dydt[mh + ke_index(k, i)] =
                  - beta_e * i * V(S, k, i, Msi) +
                  gamma * ((i + 1) * V(S, k + 1, i + 1, Msi) - i * V(S, k, i, Msi)) +
                  incidence_s_neighbors / susceptible_s_neighbors *
                      ((s + 1) * V(S, k, i - 1, Msi) - s * V(S, k, i, Msi)) -
                  internal_infection_rate / all_s * V(S, k, i, Msi) +
                  internal_infection_rate / all_s * ((s + 1) * V(S, k, i - 1, Msi) - s * V(S, k, i, Msi));
      }
  double incidence_i_neighbors = 0;
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          incidence_i_neighbors += beta_e * i * i * V(S, k, i, Msi);
      }
  double infectious_s_neighbors = 0;
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          infectious_s_neighbors += (k - i) * V(I, k, i, Msi);
      }
  // The I equations
  for (int k = 0; k <= Msi; k++)
      for (int i = 0; i <= k; i++) {
          int s = k - i;
          dydt[mh + msi + ke_index(k, i)] =
                  beta_e * i * V(S, k, i, Msi) - gamma * V(I, k, i, Msi) +
                  gamma * ((i + 1) * V(I, k + 1, i + 1, Msi) - i * V(I, k, i, Msi)) +
                  incidence_i_neighbors / infectious_s_neighbors *
                      ((s + 1) * V(I, k, i - 1, Msi) - s * V(I, k, i, Msi)) +
                  internal_infection_rate / all_s * V(S, k, i, Msi) +
                  internal_infection_rate / all_s * ((s + 1) * V(I, k, i - 1, Msi) - s * V(I, k, i, Msi));
      }
  return GSL_SUCCESS;
} </content> <use name=“equations.h”/> <use name=“ke_index.h”/>