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

#include <cmath>
 
exp(-0.25/ y_0);

double dy0dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3){
    return 8.85*pow(10, -6) - y_2*y_0*((1.2 - 1 - y_0) / (1 - y_1)) + y_0 / pow(10, -8);
}
double dy1dt(const double t, const double y_0, const double y_1, const double y_2, const double y_3){
    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){
    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){
    return pow(10, 30)*(sqrt(y_0))*(exp(-0.25 / y_0));
}
void rungeKutta_by_4(double y0, double y1, double y2, double y3, 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) >> 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);
        }
 
        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]);
        }
 
        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]);
        }
 
        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[3]);
        }
 
        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[3] + 2 * k2[3] + 2 * k3[3] + k4[3]) / 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 beg_time, max_time;
    double step;
cout << "Input begin: "; cin >> beg_time; // = 0
    cout << "Input end: "; cin >> max_time;   // = 1
cout << "Input step: "; cin >> step;      // = 0.1
rungeKutta_by_4(y_0, y_1, y_2, y_3, beg_time, max_time, step);
}