Extension SIR - 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() {

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

    	// Parameters
	double b; //Infection Rate  
	double k; //Recovery Rate

    	int numDays = 200;
	
    	std::cout << "Enter the temperature: ";
    	double temperature;
    	std::cin >> temperature;

    	if (temperature >= 0 && temperature <= 10) {
       		b = 1.0 / 2.0;
		k = 1.0 / 3.0;

    	} else if (temperature > 10 && temperature <= 20) {
        	b = 13.0 / 50.0;
		k = 13.0 / 75.0;

    	} else if (temperature > 20 && temperature <= 30) {
        	b = 2.0 / 100.0;
		k = 1.0 / 75.0;

    	} else {
        	std::cout << "Error: Enter a temperature between 0 and 30." << std::endl;

        cpgclos();
        return 1;

    	}
    
    	if (!cpgopen("/XWINDOW")) return 1;

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

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

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

    	cpgclos();

    	return 0;
}