SIR Model - C++

#include <iostream>
#include <cpgplot.h>

const double N = 7900000; // Total population
const double dt = 1.0;

double ds_dt(double s, double i, double b) {
    	return -b * s * i;
}

double di_dt(double s, double i, double b, double k) {
    	return b * s * i - k * i;
}

double dr_dt(double i, double k) {
    	return k * i;
}

void rungeKutta(double& s, double& i, double& r, double b, double k, double dt) {

    	double ds1 = ds_dt(s, i, b) * dt;
    	double di1 = di_dt(s, i, b, k) * dt;
    	double dr1 = dr_dt(i, k) * dt;

    	double ds2 = ds_dt(s + 0.5 * ds1, i + 0.5 * di1, b) * dt;
    	double di2 = di_dt(s + 0.5 * ds1, i + 0.5 * di1, b, k) * dt;
    	double dr2 = dr_dt(i + 0.5 * di1, k) * dt;

    	double ds3 = ds_dt(s + 0.5 * ds2, i + 0.5 * di2, b) * dt;
    	double di3 = di_dt(s + 0.5 * ds2, i + 0.5 * di2, b, k) * dt;
    	double dr3 = dr_dt(i + 0.5 * di2, k) * dt;

    	double ds4 = ds_dt(s + ds3, i + di3, b) * dt;
    	double di4 = di_dt(s + ds3, i + di3, b, k) * dt;
    	double dr4 = dr_dt(i + di3, k) * dt;

    	s += (ds1 + 2.0 * ds2 + 2.0 * ds3 + ds4) / 6.0;
    	i += (di1 + 2.0 * di2 + 2.0 * di3 + di4) / 6.0;
    	r += (dr1 + 2.0 * dr2 + 2.0 * dr3 + dr4) / 6.0;
}

int main() {
    
	if (!cpgopen("/XWINDOW")) return 1;

    	cpgenv(0, 140, 0, N, 0, 1);
    	cpglab("Day", "Population", "Influenza A: Hong Kong Flu");

    	// Initial conditions
    	double s = 1.0;       //Susceptible
    	double i = 1.27e-6;   //Infected
    	double r = 0.0;       //Recovered

    	// Parameters
    	double b = 1.0 / 2.0; //Contact Rate
    	double k = 1.0 / 3.0; //Recovery rate

    	// Number of days to simulate
    	int numDays = 140;

    	for (int day = 0; day < numDays; ++day) {
        	rungeKutta(s, i, r, b, k, dt);

        	cpgsci(4);
        	cpgpt1(day, N * s, 1);
		      cpgtext(10, 0.9 * N, "Susceptible");
        	cpgsci(2);
        	cpgpt1(day, N * i, 1);
		      cpgtext(10, 0.85 * N, "Infected");
        	cpgsci(3);
        	cpgpt1(day, N * r, 1);
		      cpgtext(10, 0.8 * N, "Recovered");
	}

    cpgclos();

    return 0;
}