Система ОДУ не решается

#include <iostream>
#include <functional>
#include <vector>
#include <iomanip>
 
using namespace std;
 
double dx1dt(const double t, const double x1, const double x2, const double x3) {
   return x1 - x2 + x3;
}
 
double dx2dt(const double t, const double x1, const double x2, const double x3) {
   return x1 + x2 - x3;
}
 
double dx3dt(const double t, const double x1, const double x2, const double x3) {
   return 2*x1 - 2*x2;
}
 
void rungeKutta(double x, double y, double z, double t0, double max_t, double h) {
 
   double k1[3], k2[3], k3[3], k4[3];
 
   double t = t0;
 
   vector<function<double(double, double, double, double)>> vfunc;
   vfunc.push_back(dx1dt);
   vfunc.push_back(dx2dt);
   vfunc.push_back(dx3dt);
 
   size_t max_iter = (max_t - t0) / h;
 
 
   for (int i = 0; i < max_iter; i++)
   {
      cout << t << "\t" << setprecision(5) << x << "\t" << y << "\t" << z << endl;
 
      for (int j(0); j < 3; ++j) {
         k1[j] = h*vfunc[j](t, x, y, z);
      }
 
      for (int j(0); j < 3; ++j) {
         k2[j] = h*vfunc[j](t + 0.5*h, x + 0.5*k1[0], y + 0.5*k1[1], z + 0.5*k1[2]);
      }
 
      for (int j(0); j < 3; ++j) {
         k3[j] = h*vfunc[j](t + 0.5*h, x + 0.5*k2[0], y + 0.5*k2[1], z + 0.5*k2[2]);
      }
 
      for (int j(0); j < 3; ++j) {
         k4[j] = h*vfunc[j](t + h, x + k3[0], y + k3[1], z + k3[2]);
      }
      
      x += (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.;
      y += (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.;
      z += (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.;
 
      t = t + h;
   }
   cout << t << setprecision(5) << "\t" << x << "\t" << y << "\t" << z << endl;
}
 
int main() {
   double x0 = 1;
   double y0 = 1;
   double z0 = 1;
   double beg_time = 0;
   double max_time = 1;
   double step = 0.1;
   rungeKutta(x0, y0, z0, beg_time, max_time, step);
}

double dy0dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3, const double tau, const double P, const double eps, const double alpha, const double I0){
    return eps - y_2*y_0*((P - 1 - y_0) / (1 - y_1)) + y_0 / tau;
}
double dy1dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3, const double tau, const double P, const double eps, const double alpha, const double I0){
    return y_0 * y_2;
}
double dy2dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3, const double tau, const double P, const double eps, const double alpha, const double I0){
    return y_0 * y_3;
}
 
double dy3dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3, const double tau, const double P, const double eps, const double alpha, const double I0){
    return I0*(sqrt(y_0))*pow(2.718, (-alpha / y_0)); // не видит exp = 2.718
}
 
void rungeKutta_by_4(double y0, double y1, double y2, double y3, double tau,  double P, double eps, double alpha, double I0, double t0, double max_t, double h){
    double k1[4], k2[4], k3[4], k4[4];
 
    double t = t0;
 
    vector<function<double(double, double, double, double, double, double, double, double, double, double) >> vfunc;
    vfunc.push_back(dy0dt);
    vfunc.push_back(dy1dt);
    vfunc.push_back(dy2dt);
    vfunc.push_back(dy3dt);
 
    size_t max_iter = (max_t - t0) / h;
 
    for (int i = 0; i < max_iter; i++)
    {
        cout << t << "\t" << setprecision(5) << y0 << "\t" << y1 << "\t" << y2 << "\t" << y3 << endl;
 
        for (int j(0); j < 4; ++j) {
            k1[j] = h*vfunc[j](t, y0, y1, y2, y3, tau, P, eps, alpha, I0);
        }
 
        for (int j(0); j < 4; ++j) {
            k2[j] = h*vfunc[j](t + 0.5*h, y0 + 0.5*k1[0], y1 + 0.5*k1[1], y2 + 0.5*k1[2], y3 + 0.5*k1[3], tau, P, eps, alpha, I0);
        }
 
        for (int j(0); j < 4; ++j) {
            k3[j] = h*vfunc[j](t + 0.5*h, y0 + 0.5*k2[0], y1 + 0.5*k2[1], y2 + 0.5*k2[2], y3 + 0.5*k3[3], tau, P, eps, alpha, I0);
        }
 
        for (int j(0); j < 4; ++j) {
            k4[j] = h*vfunc[j](t + h, y0 + k3[0], y1 + k3[1], y2 + k3[2], y3 + k3[4], tau, P, eps, alpha, I0);
        }
 
        y0 += (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.;
        y1 += (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.;
        y2 += (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.;
        y3 += (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.;
 
        t = t + h;
    }
    cout << t << setprecision(5) << "\t" << y0 << "\t" << y1 << "\t" << y2 << "\t" << y3 << endl;
 
}
 
int main() {
        double y_0, y_1, y_2, y_3;
    cout << "Input y0: "; cin >> y_0; // 10^(-10)
    cout << "Input y1: "; cin >> y_1; // 0
    cout << "Input y2: "; cin >> y_2; // 0
    cout << "Input y3: "; cin >> y_3; // 0
 
        // заданные константы
    double alpha;
    double I0;
    double P;
    double tau ;
 
        // зависит от t
    double eps = 8; 
 
    double beg_time, max_time;
    cout << "Input begin: "; cin >> beg_time; // = 0
    cout << "Input end: "; cin >> max_time;   // = 1
 
    double step; 
    cout << "Input step: "; cin >> step;      // = 0.1
 
    rungeKutta_by_4(y_0, y_1, y_2, y_3, beg_time, max_time, step); // сюда нужно дописывать константы 
}

//40 строка 5 параметр функции
for (int j(0); j < 4; ++j) {
   k4[j] = h*vfunc[j](t + h, y0 + k3[0], y1 + k3[1], y2 + k3[2], y3 + k3[4], tau, P, eps, alpha, I0);
}

// 50 строка
y3 += (k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.;

y3 += (k1[3] + 2 * k2[3] + 2 * k3[3] + k4[3]) / 6.;