СОДУ как сделать, если метод рунге-кутта 4го порядка не решает ее

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 x0, y0, z0, y1, y2, y3, beg_time, max_time, step;
 
    cout << "Input y0: "; cin >> y0;          // = 10^(-10) = 10e-10
    cout << "Input y1: "; cin >> y1;          // = 0
    cout << "Input y2: "; cin >> y2;          // = 0
    cout << "Input y3: "; cin >> y3;          // = 0
 
    cout << "Input begin: "; cin >> beg_time; // = 0
    cout << "Input end: "; cin >> max_time;   // = 1
 
    cout << "Input step: "; cin >> step; // = 0.1
 
    cout << "\n** \t Runge-Kutta for a system of 3 equations: \t **" << endl << endl;
    rungeKutta_by_4(y0, y1, y2, y3, beg_time, max_time, step);
 
    cout << endl;
    system("pause");
    return 0;
}

double y_0(1e-10), y_1(0), y_2(0), y_3(0);
double beg_time(0), max_time(2e-5);
double step(1e-10);

#include <iostream>
#include <iomanip>
#include <functional>
#include <vector>
#include <boost/multiprecision/cpp_bin_float.hpp>
 
using namespace boost::multiprecision;
 
typedef number<cpp_bin_float<200>> cpp_dec_float_unlim;
 
#define double cpp_dec_float_unlim
 
// для главной системы 4х уравнений
 
double dy0dt(const double& t, const double& y_0, const double& y_1, const double& y_2, const double& y_3) {
   return double(8.85e6) - y_2*y_0*((double(0.2) - y_0) / (double(1.) - y_1)) + y_0 / double(1e-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 double(1e30)*sqrt(y_0)*exp(double(-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;
 
   std::vector<std::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_t - t0) / h); i++)
   {
      std::cout << t << "\t" << std::setprecision(5) << y0 << "\t" << y1 << "\t" << y2 << "\t" << y3 << std::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 + double(0.5)*h, y0 + double(0.5)*k1[0], y1 + double(0.5)*k1[1], y2 + double(0.5)*k1[2], y3 + double(0.5)*k1[3]);
      }
 
      for (int j = 0; j < 4; ++j) {
         k3[j] = h*vfunc[j](t + double(0.5)*h, y0 + double(0.5)*k2[0], y1 + double(0.5)*k2[1], y2 + double(0.5)*k2[2], y3 + double(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] + double(2.) * k2[0] + double(2.) * k3[0] + k4[0]) / double(6.);
      y1 += (k1[1] + double(2.) * k2[1] + double(2.) * k3[1] + k4[1]) / double(6.);
      y2 += (k1[2] + double(2.) * k2[2] + double(2.) * k3[2] + k4[2]) / double(6.);
      y3 += (k1[3] + double(2.) * k2[3] + double(2.) * k3[3] + k4[3]) / double(6.);
 
      t = t + h;
   }
   std::cout << t << std::setprecision(5) << "\t" << y0 << "\t" << y1 << "\t" << y2 << "\t" << y3 << std::endl;
}
 
void foo() {
   std::cout << "\nfoo\n";
 
   double y_0(1e-10), y_1(0), y_2(0), y_3(0);
 
   std::cout << y_0;
 
   double beg_time(0), max_time(2e-5);
   double step(1e-11);
 
   rungeKutta_by_4(y_0, y_1, y_2, y_3, 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) {
   std::cout << double(8.85e6) - y_2*y_0*((double(0.2) - y_0) / (double(1.) - y_1)) + y_0 / double(1e-8);
   return double(8.85e6) - y_2*y_0*((double(0.2) - y_0) / (double(1.) - y_1)) + y_0 / double(1e-8);
}