Go back to Flu Mortality

#include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_odeiv.h>

int func (double t, const double y[], double dydt[], void *params);

void SIR_solver (double *params, double *t, double *yin, double *ansS) {

const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);

gsl_odeiv_system sys = {func, NULL, 2, params};

double t1 = t[0];
double t2 = t[1];

double h = 1e-6;

double y1 = yin[0];
double y2 = yin[1];

double y[2] = { y1, y2 };
  
  for (int i = 1; i <= t2; i++)
       {
         double ti = i * t1 / t2;
   
         while (t1 < ti)
           {
             gsl_odeiv_evolve_apply (e, c, s, &sys, &t1, ti, &h, y);
           }
           
         ansS[i] = y[0];
       }

gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);

}

int func (double t, const double y[], double dydt[], void *params) {

double beta = ((double *)params)[0];
double alpha = ((double *)params)[1];
double N = ((double *)params)[2];

dydt[0] = -1*beta*y[0]*y[1];
dydt[1] = (beta*y[0]*y[1]) - (alpha*y[1]);

return GSL_SUCCESS;

}