Program stops after random number of iterations

I'm making a program which allow me to find the numerical solution of a 2-body problem involving the Sun and the Earth by using a Runge-Kutta scheme. I'm using a structure to define the solution properties, another one to define the numerical method used and a class to define the planet properties. The result is finally wrote on a file. I have no problem in running the program with every Runge-Kutta method defined but the Runge-Kutta of order 7 (RK7). The problem is in the cycle inside the RK function of the main file: after a random number of iterations (usually no more than i = 4) the program stops ("Main.exe stopped working"). I checked what could be the reason of that and it seems the problem is at

supp = f(fa, fb, Ms);

which is a support vector that receives the result of the 2-body system of differential equations (function F2 for instance). I don't understand what could be the reason since the program works fine in any other case (I already checked the RK7 Butcher tableau values). I don't know if this can help, but trying to print the values of the p vector from the F2 function, I noticed that, after the aforementioned random number of iterations, the program stops sometimes without printing all the values of p.

Main.cpp

#include<iostream>
#include<math.h>
#include<vector>
#include "celestial_body.cpp"
#include "solution.cpp"
#include "RK_schemes.cpp"

using namespace std;

const double P = 365.25;
const double G = 1.000009448418457; 
const double Ms = 1.0;
const double pi = 3.14159265358979323846;

typedef vector<double> (* fun) (double, vector<double>, double);

vector<double> F1(double t, vector<double> x, double m) {
    vector<double> p(4,0);
    //r system
    p[0] = x[1];
    p[1] = x[0]*pow(x[3],2) - (G*m)/pow(x[0],2);
    //theta system
    p[2] = x[3];
    p[3] = -2*x[1]*x[3]/x[0];
    cout << p[0] << " " << p[1] << " " << p[2] << " " << p[3] << endl;
    return p;
}

vector<double> F2(double t, vector<double> x, double m) {
    vector<double> p(4,0);
    double d = sqrt(pow(x[0],2) + pow(x[2],2));

    //x system
    p[0] = x[1];
    p[1] = -(G*m*(x[0]))/pow(d,3);
    //y system
    p[2] = x[3];
    p[3] = -(G*m*x[2])/pow(d,3);

    return p;
}

void RK(fun f, CB ob1, int N, double tf, double t0 = 0) {
    double fa, h = (tf - t0)/N;
    cout << "h is: " << h << endl;

    Sol U(N, 4); //U.dim = 4
    U.set_starting_values(t0, ob1.CB::getparam()); 

    //RK_model RK(4);
    RK_model RK(7);
    //RK.set_RK4();
    RK.set_RK7();

    vector<double> sum1(U.dim,0);
    vector<double> sum2(U.dim,0);
    vector<double> K(U.dim*U.dim,0);
    vector<double> fb(U.dim,0);
    vector<double> supp(U.dim,0);

    int k = 0;
    for (int n = 0; n < N-1; n++) {
        U.t[n+1] = t0 + (n+1)*h;
        for (int i = 0; i < RK.nu; i++) {
            fa = U.t[n] + RK.c[i]*h;
            for (int j = 0; j < RK.nu; j++) {
                if (i == j || i < j) {
                    break;
                }
                while (k < U.dim) {
                    sum1[k] += RK.a[j*RK.nu + i]*K[k*RK.nu + j];
                    k++;
                }
                k = 0;
            }
            while (k < U.dim) {
                fb[k] = U.u[k*N+n] + h*sum1[k];
                k++;
            }
            k = 0;

            fill(sum1.begin(), sum1.end(), 0);
            supp = f(fa, fb, Ms);

            while (k < U.dim) {
                K[k*RK.nu + i] = supp[k];
                sum2[k] += RK.b[i]*K[k*RK.nu + i];
                k++;
            }
            k = 0;
        }
        while (k < U.dim) {
            U.u[k*N+n+1] = U.u[k*N+n] + h*sum2[k];
            k++;
        }
        k = 0;

        fill(sum2.begin(), sum2.end(), 0);

    }
    
    U.write_to_file();
    U.del();
    RK.del();
    vector<double>().swap(sum1);
    vector<double>().swap(sum2);
    vector<double>().swap(K);
    vector<double>().swap(fb);

}

int main() {
    CB E(3.002513826043238e-06);
    int fn;
    double N, tf;
        
    cout << "Choose function (1: polar; 2: cartesian): ";
    cin >> fn;
    assert(fn == 1 || fn == 2);
    cout << "Choose final time (starting time have been set to zero): ";
    cin >> tf;
    cout << "Choose number of iterations: ";
    cin >> N;

    //Set starting position and velocity of Earth
    if (fn == 1) {
        //start apocenter
        //E.CB::setpos(1.016719016709908, 0);
        //E.CB::setvel(0, 0.016638215898855);

        //start pericenter
        E.CB::setpos(0.983249288988710, 0);
        E.CB::setvel(0, 0.017791966932260);

        RK(F1, E, N, tf);
    } else if(fn == 2) {
        //start apocenter
        //E.CB::setpos(1.016719016709908, 0);
        //E.CB::setvel(0, 0.983372498367391);

        //start pericenter
        E.CB::setpos(0.983249288988710, 0);
        E.CB::setvel(0, 1.016789362554682);

        RK(F2, E, N, tf);
    }

solution.cpp

#ifndef SOLUTION_4D
#define SOLUTION_4D

#include<iostream>
#include<fstream>
#include<vector>

using namespace std;

struct Sol {
    vector<double> u;
    vector<double> t;
    int N;
    int dim;

    Sol(int n, int d) {
        dim = d;
        N = n;
        for (int i = 0; i < dim*N; i++) {
            u.push_back(0);
            if (i < N) {
                t.push_back(0);
            }
        }
    }

    void set_starting_values(double t0, vector<double> p) {
        t[0] = t0;
        for (int i = 0; i < dim; i++) {
            u[i*N] = p[i];
        }
    } 

    void write_to_file() {
        ofstream myfile("result.m", ios_base::out);
        myfile.setf(ios_base::left, ios_base::adjustfield);
        myfile.setf(ios_base::scientific);
        myfile.precision(15); 
        int k = 0;
        myfile << "u" << " = [ " << endl;
        for (int i = 0; i < N; i++) {
            myfile << t[i] << " ";
            while (k < dim) {
                if (k == dim-1) {
                    myfile << u[k*N+i] << endl;
                } else {
                    myfile << u[k*N+i] << " ";
                }
                k++;
            }
            k = 0;
        }

        myfile << "]; \n";
        myfile.close();

    }

    void del() {
        vector<double>().swap(t);
        vector<double>().swap(u);;
    }

};

#endif

celestial_body.cpp

#ifndef CELESTIAL_BODY
#define CELESTIAL_BODY

#include<vector>

using namespace std;

class CB {
private:
    double M;
    vector<double> param; //contain xpos, xvel, ypos, yvel in this order

public:
    CB(double m) {
        M = m;
        for (int i = 0; i < 4; i++) {
            param.push_back(0);
        }
    }

    void setpos(double x, double y) {
        param[0] = x;
        param[2] = y;
        }

    vector<double> getpos() {
        vector<double> pos(2,0);
        pos[0] = param[0];
        pos[1] = param[2];
        return pos;
    }

    void setvel(double vx, double vy) {
        param[1] = vx;
        param[3] = vy;
    }

    vector<double> getvel() {
        vector<double> vel(2,0);
        vel[0] = param[1];
        vel[1] = param[3];
        return vel;
    }

    vector<double> getparam() {
        return param;
    }

    double getmass() {
        return M;
    }

};

#endif

RK_schemes.cpp

#ifndef RK_SCHEMES
#define RK_SCHEMES

#include<cassert>
#include<iomanip>
#include<vector>
#include<math.h>


struct RK_model {
    vector<double> c;
    vector<double> b;
    vector<double> a;
    int nu;

    RK_model(int stages) {
        nu = stages;

        for (int i = 0; i < nu*nu; i++) {
            if (i < nu) {
                c.push_back(0);
                b.push_back(0);
            }
            a.push_back(0);
        } 
    } 

    void set_EE() {
        assert(nu == 1);
        b[0] = 1;
    }

    void set_Heun() {
        assert(nu == 2);
        c[1] = 1.0;

        b[0] = 1.0/2.0;
        b[1] = 1.0/2.0;

        a[1] = 1.0;
    }

    void set_modified_Euler() {
        assert(nu == 2);
        c[1] = 1.0/2.0;

        b[1] = 1.0;

        a[1] = 1.0/2.0;
    }

    void set_Shu() {
        assert(nu == 3);
        c[1] = 1.0/3.0;
        c[2] = 2.0/3.0;

        b[0] = 1.0/4.0;
        b[2] = 3.0/4.0;

        a[1] = 1.0/3.0;
        a[5] = 2.0/3.0;
    }

    void set_RK4() {
        assert(nu == 4);
        c[1] = 1.0/2.0;
        c[2] = 1.0/2.0;
        c[3] = 1.0;

        b[0] = 1.0/6.0;
        b[1] = 1.0/3.0;
        b[2] = 1.0/3.0;
        b[3] = 1.0/6.0;

        a[1] = 1.0/2.0;
        a[6] = 1.0/2.0;
        a[11] = 1.0;
    }

    void set_RK7() {
        assert(nu == 7);
        double r = sqrt(21);
        c[1] = 1.0;
        c[2] = 0.5;
        c[3] = 2.0/3.0;
        c[4] = (7.0 - r)/14.0;
        c[5] = (7.0 + r)/14.0;
        c[6] = 1.0;

        b[0] = 1.0/20.0;
        b[2] = 16.0/45.0;
        b[4] = 49.0/180.0;
        b[5] = 49.0/180.0;
        b[6] = 1.0/20.0;

        a[1] = 1;
        a[2] = 3.0/8.0;
        a[3] = 8.0/27.0;
        a[4] = 3*(3*r - 7)/392.0;
        a[5] = -(231 + 51*r)/392.0;
        a[6] = (22 + 7*r)/12.0;

        a[9] = 1.0/8.0;
        a[10] = 2.0/27.0;
        a[11] = (7 - r)/49.0;
        a[12] = -(7 + r)/49.0;
        a[13] = 2.0/3.0;

        a[17] = 8.0/27.0;
        a[18] = 6*(7 - r)/49.0;
        a[19] = -40*r/245.0;
        a[20] = 2*(7*r - 5)/9.0;

        a[25] = 3*(21 - r)/392.0;
        a[26] = 3*(21 + 121*r)/1960.0;
        a[27] = -63*(3*r - 2)/180.0;

        a[33] = (6 + r)/5.0;
        a[34] = -7*(49 + 9*r)/90.0;

        a[41] = 7*(7 - r)/18.0;
    }

    void del() {
        vector<double>().swap(a);
        vector<double>().swap(b);
        vector<double>().swap(c);
    }

    void print_tab() {
        cout << endl;
        cout << "Butcher tableau" << endl;
        for (int i = 0; i < nu; i++) {
            cout << c[i] << setprecision(2) << " " << "|" << " ";
            for (int j = 0; j < nu; j++) {
                cout << setw(3) << left << a[j*nu + i];
            }
            cout << endl;
        }

        for (int i = 0; i < nu + 5; i++) {
            cout << "-";
        }
        cout << endl;
        for (int i = 0; i < nu; i++) {
            cout << b[i] << " ";
        }
        cout << endl;
    }
};

#endif