C++ : vector subscript out of range?

I'm new to C++ and I'm trying to write a program which opens a file "parameters.txt", which has 8 space-delimited numbers in it. Then I use these numbers in a separate function RK4() to get 3 arrays (y1, y2, t). Then I manipulate these arrays in main() and output a file "output.txt" to the working directory

This is my code:

#include <string>
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
#include <algorithm>
#include <iterator>
#include <tuple>
using namespace std;

tuple<vector<double>, vector<double>, vector<double>> RK4() {

    //open parameters.txt, put data into a vector
    ifstream fin("parameters.txt");
    vector<double> data;
    int element;
    while (fin >> element) {
        data.push_back(element);
    }

    //define tspan
    vector<double> tspan(2);
    tspan[0] = 0.0;
    tspan[1] = data[7];

    //define y0
    vector<double> y0(4);
    y0[0] = data[4];
    y0[1] = data[5];
    y0[2] = 0.0;
    y0[3] = 0.0;
    double theta1 = y0[0];
    double theta2 = y0[1];
    double omega1 = y0[2];
    double omega2 = y0[3];

    //define stepSize
    double stepSize;
    stepSize = data[6];

    //define range
    int range = int(tspan[1] / stepSize);

    //define other constants
    double m1, m2, l1, l2;
    m1 = data[0];
    m2 = data[1];
    l1 = data[2];
    l2 = data[3];
    double g = 9.81;

    //define y, t vectors
    vector<double> y1(range);
    vector<double> y2(range);
    vector<double> y3(range);
    vector<double> y4(range);
    vector<double> t(range);
    for (double i = 0.0; i < 1.0 * range; i++) {
        t[i] = i * stepSize;
    }

    //enter y0 into first value
    y1[0] = theta1;
    y2[0] = theta2;
    y3[0] = omega1;
    y4[0] = omega2;

    //loop to find y, t vectors
    for (int i = 0; i < range - 1; i++) {
        //finding all k values:
        //k1
        double dTheta1_1 = y3[i];
        double dOmega1_1 = (-g * (2 * m1 + m2) * sin(y1[i]) - m2 * g * sin(y1[i] - 2 * y2[i]) - 2 * sin(y1[i] - y2[i]) * m2 * (pow(y4[i], 2) * l2 + pow(y3[i], 2) * l1 * cos(y1[i] - y2[i]))) / (l1 * (2 * m1 + m2 - m2 * cos(2 * y1[i] - 2 * y2[i])));
        double dTheta2_1 = y4[i];
        double dOmega2_1 = (2 * sin(y1[i] - y2[i]) * (pow(y3[i], 2) * l1 * (m1 + m2) + g * (m1 + m2) * cos(y1[i]) + pow(y4[i], 2) * l2 * m2 * cos(y1[i] - y2[i]))) / (l2 * (2 * m1 + m2 - m2 * cos(2 * y1[i] - 2 * y2[i])));

        //k2
        double dTheta1_2 = y3[i] + 0.5 * stepSize * dTheta1_1;
        double dOmega1_2 = (-g * (2 * m1 + m2) * sin(y1[i] + 0.5 * stepSize * dTheta1_1) - m2 * g * sin((y1[i] + 0.5 * stepSize * dTheta1_1) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_1)) - 2 * sin((y1[i] + 0.5 * stepSize * dTheta1_1) - (y2[i] + 0.5 * stepSize * dTheta2_1)) * m2 * (pow(y4[i] + 0.5 * stepSize * dOmega2_1, 2) * l2 + pow(y3[i] + 0.5 * stepSize * dOmega1_1, 2) * l1 * cos((y1[i] + 0.5 * stepSize * dTheta1_1) - (y2[i] + 0.5 * stepSize * dTheta2_1)))) / (l1 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + 0.5 * stepSize * dTheta1_1) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_1))));
        double dTheta2_2 = y4[i] + 0.5 * stepSize * dTheta2_1;
        double dOmega2_2 = (2 * sin((y1[i] + 0.5 * stepSize * dTheta1_1) - (y2[i] + 0.5 * stepSize * dTheta2_1)) * (pow(y3[i] + 0.5 * stepSize * dOmega1_1, 2) * l1 * (m1 + m2) + g * (m1 + m2) * cos(y1[i] + 0.5 * stepSize * dTheta1_1) + pow(y4[i] + 0.5 * stepSize * dOmega2_1, 2) * l2 * m2 * cos((y1[i] + 0.5 * stepSize * dTheta1_1) - (y2[i] + 0.5 * stepSize * dTheta2_1)))) / (l2 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + 0.5 * stepSize * dTheta1_1) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_1))));

        //k3
        double dTheta1_3 = y3[i] + 0.5 * stepSize * dTheta1_2;
        double dOmega1_3 = (-g * (2 * m1 + m2) * sin(y1[i] + 0.5 * stepSize * dTheta1_2) - m2 * g * sin((y1[i] + 0.5 * stepSize * dTheta1_2) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_2)) - 2 * sin((y1[i] + 0.5 * stepSize * dTheta1_2) - (y2[i] + 0.5 * stepSize * dTheta2_2)) * m2 * (pow(y4[i] + 0.5 * stepSize * dOmega2_2, 2) * l2 + pow(y3[i] + 0.5 * stepSize * dOmega1_2, 2) * l1 * cos((y1[i] + 0.5 * stepSize * dTheta1_2) - (y2[i] + 0.5 * stepSize * dTheta2_2)))) / (l1 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + 0.5 * stepSize * dTheta1_2) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_2))));
        double dTheta2_3 = y4[i] + 0.5 * stepSize * dTheta2_2;
        double dOmega2_3 = (2 * sin((y1[i] + 0.5 * stepSize * dTheta1_2) - (y2[i] + 0.5 * stepSize * dTheta2_2)) * (pow(y3[i] + 0.5 * stepSize * dOmega1_2, 2) * l1 * (m1 + m2) + g * (m1 + m2) * cos(y1[i] + 0.5 * stepSize * dTheta1_2) + pow(y4[i] + 0.5 * stepSize * dOmega2_2, 2) * l2 * m2 * cos((y1[i] + 0.5 * stepSize * dTheta1_2) - (y2[i] + 0.5 * stepSize * dTheta2_2)))) / (l2 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + 0.5 * stepSize * dTheta1_2) - 2 * (y2[i] + 0.5 * stepSize * dTheta2_2))));

        //k4
        double dTheta1_4 = y3[i] + stepSize * dTheta1_3;
        double dOmega1_4 = (-g * (2 * m1 + m2) * sin(y1[i] + stepSize * dTheta1_3) - m2 * g * sin((y1[i] + stepSize * dTheta1_3) - 2 * (y2[i] + stepSize * dTheta2_3)) - 2 * sin((y1[i] + stepSize * dTheta1_3) - (y2[i] + stepSize * dTheta2_3)) * m2 * (pow(y4[i] + stepSize * dOmega2_3, 2) * l2 + pow(y3[i] + stepSize * dOmega1_3, 2) * l1 * cos((y1[i] + stepSize * dTheta1_3) - (y2[i] + stepSize * dTheta2_3)))) / (l1 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + stepSize * dTheta1_3) - 2 * (y2[i] + stepSize * dTheta2_3))));
        double dTheta2_4 = y4[i] + stepSize * dTheta2_3;
        double dOmega2_4 = (2 * sin((y1[i] + stepSize * dTheta1_3) - (y2[i] + stepSize * dTheta2_3)) * (pow(y3[i] + stepSize * dOmega1_3, 2) * l1 * (m1 + m2) + g * (m1 + m2) * cos(y1[i] + stepSize * dTheta1_3) + pow(y4[i] + stepSize * dOmega2_3, 2) * l2 * m2 * cos((y1[i] + stepSize * dTheta1_3) - (y2[i] + stepSize * dTheta2_3)))) / (l2 * (2 * m1 + m2 - m2 * cos(2 * (y1[i] + stepSize * dTheta1_3) - 2 * (y2[i] + stepSize * dTheta2_3))));


        double theta1New = y1[i] + (stepSize / 6.0) * (dTheta1_1 + 2 * dTheta1_2 + 2 * dTheta1_3 + dTheta1_4);
        double omega1New = y3[i] + (stepSize / 6.0) * (dOmega1_1 + 2 * dOmega1_2 + 2 * dOmega1_3 + dOmega1_4);
        double theta2New = y2[i] + (stepSize / 6.0) * (dTheta2_1 + 2 * dTheta2_2 + 2 * dTheta2_3 + dTheta2_4);
        double omega2New = y4[i] + (stepSize / 6.0) * (dOmega2_1 + 2 * dOmega2_2 + 2 * dOmega2_3 + dOmega2_4);


        // updating y arrays 
        y1[i + 1] = theta1New;
        y2[i + 1] = theta2New;
        y3[i + 1] = omega1New;
        y4[i + 1] = omega2New;

    }
    return make_tuple(y1, y2, t);
}

int main() {

    //open parameters.txt, put data into a vector
    ifstream fin("parameters.txt");
    vector<double> data;
    int element;
    while (fin >> element) {
        data.push_back(element);
    }
   
    //define tspan
    vector<double> tspan(2);
    tspan[0] = 0.0;
    tspan[1] = 10.0;

    //define stepSize
    double stepSize;
    stepSize = data[6];

    //define range
    int const range = 1000;

    //define other constants
    double l1 = data[2];
    double l2 = data[3];

    //get y1, y2, t from RK4 function
    auto temp = RK4();
    vector<double> y1 = get<0>(temp);
    vector<double> y2 = get<1>(temp);
    vector<double> t = get<2>(temp);
    double x_1[range], y_1[range], x_2[range], y_2[range];

    //define x_1, x_2, y_1, y_2
    for (int i = 0; i < range; i++) {
        x_1[i] = { sin(y1[i]) * l1 };
        y_1[i] = { -cos(y1[i]) * l1 };
        x_2[i] = { sin(y1[i]) * l1 + sin(y2[i]) * l2 };
        y_2[i] = { -cos(y1[i]) * l1 - cos(y2[i]) * l2 };
    }

    //writing x,y positions at time t to output.txt
    ofstream myfile;
    myfile.open("C:\\mydirectory\\output.txt");
    if (myfile.is_open()) {
        myfile << "t: " << endl;
        for (int i = 0; i < range; i++) {
            myfile << t[i] << " ";
        }
        cout << endl;
        myfile << "x_1: " << endl;
        for (int i = 0; i < range; i++) {
            myfile << x_1[i] << " ";
        }
        cout << endl;
        myfile << "y_1: " << endl;
        for (int i = 0; i < range; i++) {
            myfile << y_1[i] << " ";
        }
        cout << endl;
        myfile << "x_2: " << endl;
        for (int i = 0; i < range; i++) {
            myfile << x_2[i] << " ";
        }
        cout << endl;
        myfile << "y_2: " << endl;
        for (int i = 0; i < range; i++) {
            myfile << y_2[i] << " ";
        }
        cout << endl;
        myfile.close();

    }
    else cout << "Unable to open file";

    return 0;
}

When I try to build and run the program (in Visual Studio) I get this error:

The thread 0x22c0 has exited with code 0 (0x0). Debug Assertion Failed!

Expression: vector subscript out of range

For information on how your program can cause an assertion failure, see the Visual C++ documentation on asserts.

The program has exited with code 3 (0x3).

When I try to debug the program, I get:

Exception Thrown c++ coursework new.exe has triggered a breakpoint.

What is the problem?