Вылетает программа

void func(double y0, double y1, double y2, double y3, double t0, double max_t, double h, int n, double matr[N][N]) {
 
    double t = t0;
    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;
        y0 += matr_inv[N][N]*h*dy0dt(y0,y1,y2,y3,t0);
        y1 += (matr_inv[N][N] * h*dy1dt(y0, y1, y2, y3, t0));
        y2 += (matr_inv[N][N] * h*dy2dt(y0, y1, y2, y3, t0));
        y3 += (matr_inv[N][N] * h*dy3dt(y0, y1, y2, y3, t0));
        t = t + h;
    }
    std::cout << t << "\t" << std::setprecision(5) << y0 << "\t" << y1 << "\t" << y2 << "\t" << y3 << std::endl;
}

void dydt(double y[N],double f[N], double t) 
{
    f[0]= 8.85e4 - y[2] * y[0]*(0.2 - y[0]) / (1. - y1))+ y[0] / 1e-8;
    f[1]=y[0] * y[2];
    f[2]=y[0] * y[3];
    f[3]=1e30*sqrt(y[0])*exp(-0.25) / y0;
}

void mul_1(double matr[N][N], double v[N], double result[N]) {
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            result[i] += matr[i][j] * v[j];
        }
    }
}
 
void mul(double v[N], double h, double result[N]) {
    for (int i = 0; i < N; i++) {
            v[i] *= h;
            result[i] = v[i];
    }
}
 
void add(double v1[N], double v2[N], double result[N]) {
    for (int i = 0; i < N; i++) {
            v1[i] += v2[i];
            result[i] = v1[i];
    }
}

void mul_1(double matr[N][N], double v[N], double result[N]) {
    for (int i = 0; i < N; ++i) {
        rsult[i]=0;
        for (int j = 0; j < N; ++j) {
            result[i] += matr[i][j] * v[j];
        }
    }
}
void mul(double v[N], double h, double result[N]) {
    for (int i = 0; i < N; i++) {
             result[i] = v[i]*h;
    }
}
void add(double v1[N], double v2[N], double result[N]) {
    for (int i = 0; i < N; i++) {
            result[i] = v1[i]+v2[i];
    }
}

void rosenbrock(double y[N], double t0, double max_t, double h, int n, double matr_inv[N][N]) {
     double t = t0;
     double ynew[N];
     double f[N];
     for t = t0; t <= max_t ;t+=h)
     {
       std::cout << t << "\t" << std::setprecision(5) << y[0] << "\t" << y[1] << "\t" << y[2] << "\t" << y[3] << std::endl;
       dydt(y,f, t);
       mul(matr_inv,f,ynew);
       mul(ynew,h,ynew);
       add(ynew,y,y);
     }
}

void ros(double y[N], double t0, double max_t, double h) {
    double t = t0;
    double ynew[N];
    double f[N];
    double adj[N][N];  // To adjoint of A[][] 
    double inv[N][N]; // To store inverse of A[][] 
 
    for (t = t0; t <= max_t; t += h) {
        
        //единичная матрица
        double a[N][N];
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j) a[i][j] = 1;
                else a[i][j] = 0;
            }
        }
        
        // якобиан системы
        double jac[N][N] = {
                            {(2 * (y[2] * (y[0] - 0.1) - 50000000.0 * (y[1] - 1.0))) / (1.0 - y[1]), ((y[0] - 0.2)*y[0] * y[2]) / pow((1.0 - y[1]), 2.), ((y[0] - 0.2)*y[0]) / (1.0 - y[1]),  0.0 },
                            {y[2], 0.0, y[0], 0.0},
                            {y[3], 0.0, 0.0, y[0]},
                            {(5 * pow(10, 29.)*(exp(-0.75 / y[0]))*(1.5 + y[0])) / pow(y[0], 3. / 2.), 0.0, 0.0, 0.0}
        };
 
        //умножение матрицы на число
        double umn[N][N];
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                jac[i][j] *= (h / 2);
                umn[i][j] = jac[i][j];
            }
        }
 
        // вычитание матрицы
        double temp[N][N];
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                temp[i][j] = a[i][j] - umn[i][j];
            }
        }
 
        adjoint(temp, adj);
 
        if (inverse(temp, inv))
            PrintMatr(inv, N);
 
        std::cout << t << "\t" << std::setprecision(5) << y[0] << "\t" << y[1] << "\t" << y[2] << "\t" << y[3] << std::endl;
        dydt(y, f, t);
        mul_1(inv, f, ynew);
        mul_2(ynew, h, ynew);
        add(ynew, y, y);
    }
    std::cout << t << "\t" << std::setprecision(5) << y[0] << "\t" << y[1] << "\t" << y[2] << "\t" << y[3] << std::endl;
}
 
int main()
{
    double y[N], f[N];
    y[0] = 10e-10, y[1] = 0, y[2] = 0, y[3] = 0;
    double h = 1e-9;
    double beg_time = 0.0, max_time = 2.5;
 
    cout << "Ros method:\n";
    ros(y, beg_time, max_time, h);
 
    system("pause");
 
    return 0;
}

SetJak(double y[N],double jac[N][N])
{
    jac[0][0]=(2. * (y[2] * (y[0] - 0.1) - 5.e7 * (y[1] - 1.0))) / (1.0 - y[1]);
    jac[0][1]=((y[0] - 0.2)*y[0] * y[2]) / pow((1.0 - y[1]), 2.);
    jak[0][2]= ((y[0] - 0.2)*y[0]) / (1.0 - y[1]);
    jak[0][3]=0.;
    jak[1][0]=y[2];
    jak[1][1]=0.;
    jak[1][2]=y[0];
    jak[1][3]=0.;
    jak[2][0]=y[3];
    jak[2][1]=0.;
    jak[2][2]=0.;
    jak[2][3]=y[0];
 
    jak[3][0]=5.e29*(exp(-0.75 / y[0]))*(1.5 + y[0])) / pow(y[0], 3. / 2.);
    jak[3][1]=0.;
    jak[3][2]=0.;
    jak[3][3]=0.;
 }