Проблема с выводом в консоль выдаёт точку остановы c кодом ошибки c0000374

#include <iostream> //Подключаем все необходимые библиотеки
#include <conio.h>
#include <math.h>
#include <cmath>
#include <Windows.h>
#include <iomanip> 
using namespace std;
#define eps 0.001
 
double dU(double y)
{
    return y; //Делаем первую замену переменных вида (u'= y)
}
 
double dY(double z)
{
    return z; //Делаем вторю замену переменных вида (y'= z)
}
 
double dZ(double x, double u)
{
    return  exp(x) + u; //Делаем третию замену переменных и переносим оставшиеся слагаемые на вправо в итоге получаем (z'= e^x -xu)
}
 
double func(double v)
{
    return ((sin((sqrt(3) * v) / 2)) / (3 * sqrt(3)) + cos((sqrt(3) * v) / 2)) / (exp(v / 2)) + (v * exp(v)) / 3; //Функция для проверки вычислений
}
 
int main()
{
    SetConsoleCP(1251); //Ставим русскую локализацию на ввод данных
    SetConsoleOutputCP(1251); //Ставим русскую локализацию на вывод данных
 
    double n = 6.0; //Инициализируем константу, которая отвечает за число разбиения отрезка от 0 до 2, её можно менять
 
    double* X = new double[n + 1], *U = new double[n + 1], *Y = new double[n + 1], *Z = new double[n + 1];
    X[0] = 0, U[0] = 1, Y[0] = 0, Z[0] = 0; //Определяем начальные условия
 
    double* X1 = new double[2 * n + 2], * U1 = new double[2 * n + 2], * Y1 = new double[2 * n + 2], * Z1 = new double[2 * n + 2];
    X1[0] = 0, U1[0] = 1, Y1[0] = 0, Z1[0] = 0; //Определяем начальные условия
 
    double Xmin = 0; //Левая граница знаений X
    double Xmax = 2; //Правая граница значений X
    double hX = 0;  //Шаг с которым смещаемся по оси X
    double hX1 = 0;  //Шаг с которым смещаемся по оси X
 
    double k1U = 0, k2U = 0, k3U = 0, k4U = 0; //Объявляем все необходимые переменные 
    double k1Y = 0, k2Y = 0, k3Y = 0, k4Y = 0;
    double k1Z = 0, k2Z = 0, k3Z = 0, k4Z = 0;
 
    for (int i = Xmin; i <= n; i++) {
        U[i] = 1;
        //cout <<"U["<<i<<"]" << U[i] << endl;
    }
    for (int i = Xmin; i <= 2*n; i++) {
        U1[i + 1] = 2;
        //cout << "U1[" << i << "]" << U1[i] << endl;
    }
 
    cout << "Значение X" << "\t\t" << "Значение U(X)" << "\t\t"<< "Значение F(x)" << endl; //Выводим заголовок для последующей таблицы значений
    for (int k = n; k >= 0; k--) {
        if (abs(U1[k] - U[k/2]) > eps) {
            hX = (Xmax - Xmin) / n;  //Шаг с которым смещаемся по оси X
            for (int i = Xmin; i <= n; i++){
                //Далее Реализуем метод Рунге-Котте 4 порядка точности, только для 3 порядка дифференциального уравнения
                X[i] = (Xmin + hX * i);//Увеличиваем X c шагом hX начиная с левого конца
                k1Z = dZ(X[i], U[i]); //Считаем 1 коэфициент для Z
                k1Y = dY(Z[i]); //Считаем 1 коэфициент для Y
                k1U = dU(Y[i]); //Считаем 1 коэфициент для U
 
                k2Z = dZ(X[i] + hX / 2, U[i] + (k1U * hX) / 2); //Считаем 1 коэфициент для Z
                k2Y = dY(Z[i] + (k1Z * hX) / 2); //Считаем 2 коэфициент для Y
                k2U = dU(Y[i] + (k1Y * hX) / 2); //Считаем 2 коэфициент для U
 
                k3Z = dZ(X[i] + hX / 2, U[i] + (k2U * hX) / 2); //Считаем 3 коэфициент для Z
                k3Y = dY(Z[i] + (k2Z * hX) / 2); //Считаем 3 коэфициент для Y
                k3U = dU(Y[i] + (k2Y * hX) / 2); //Считаем 3 коэфициент для U
 
                k4Z = dZ(X[i] + hX, U[i] + k3U * hX); //Считаем 4 коэфициент для Z
                k4Y = dY(Z[i] + hX * k3Z); //Считаем 4 коэфициент для Y
                k4U = dU(Y[i] + hX * k3Y); //Считаем 4 коэфициент для U
 
                Z[i + 1] = Z[i] + (k1Z + 2 * k2Z + 2 * k3Z + k4Z) * (hX / 6); //Считаем следующее значение Z (изначально он считается от начального Z)
                Y[i + 1] = Y[i] + (k1Y + 2 * k2Y + 2 * k3Y + k4Y) * (hX / 6); //Считаем следующее значение Y (изначально он считается от начального Y)
                U[i + 1] = U[i] + (k1U + 2 * k2U + 2 * k3U + k4U) * (hX / 6); //Считаем следующее значение U (изначально он считается от начального U)
            }
            n *= 2;
            hX1 = (Xmax - Xmin) / n;
            for (int i = Xmin; i <= n; i++){
                //Далее Реализуем метод Рунге-Котте 4 порядка точности, только для 3 порядка дифференциального уравнения
                X1[i] = (Xmin + hX1 * i);//Увеличиваем X c шагом hX начиная с левого конца
                k1Z = dZ(X1[i], U1[i]); //Считаем 1 коэфициент для Z
                k1Y = dY(Z1[i]); //Считаем 1 коэфициент для Y
                k1U = dU(Y1[i]); //Считаем 1 коэфициент для U
 
                k2Z = dZ(X1[i] + hX1 / 2, U1[i] + (k1U * hX1) / 2); //Считаем 1 коэфициент для Z
                k2Y = dY(Z1[i] + (k1Z * hX1) / 2); //Считаем 2 коэфициент для Y
                k2U = dU(Y1[i] + (k1Y * hX1) / 2); //Считаем 2 коэфициент для U
 
                k3Z = dZ(X1[i] + hX1 / 2, U1[i] + (k2U * hX1) / 2); //Считаем 3 коэфициент для Z
                k3Y = dY(Z1[i] + (k2Z * hX1) / 2); //Считаем 3 коэфициент для Y
                k3U = dU(Y1[i] + (k2Y * hX1) / 2); //Считаем 3 коэфициент для U
 
                k4Z = dZ(X1[i] + hX1, U1[i] + k3U * hX1); //Считаем 4 коэфициент для Z
                k4Y = dY(Z1[i] + hX1 * k3Z); //Считаем 4 коэфициент для Y
                k4U = dU(Y1[i] + hX1 * k3Y); //Считаем 4 коэфициент для U
 
                Z1[i + 1] = Z1[i] + (k1Z + 2 * k2Z + 2 * k3Z + k4Z) * (hX1 / 6); //Считаем следующее значение Z (изначально он считается от начального Z)
                Y1[i + 1] = Y1[i] + (k1Y + 2 * k2Y + 2 * k3Y + k4Y) * (hX1 / 6); //Считаем следующее значение Y (изначально он считается от начального Y)
                U1[i + 1] = U1[i] + (k1U + 2 * k2U + 2 * k3U + k4U) * (hX1 / 6); //Считаем следующее значение U (изначально он считается от начального U)
            }
        }
        else {
            for (int k = 0; k <= n; k++){
                cout << fixed << setprecision(3) << "X[" << k << "] = " << X1[k] << "\t\t" << setprecision(5) << "U[" << k << "] = " << U1[k] << "\t\t" << "func[" << k << "] = " << func(X1[k]) << endl; //Выводим значения X и U
            }
            return 0;
        }
        k = n+1;
    }
 
    return 0;
}

//double* X = new double[n + 1]
std::vector<double> X(n + 1)
//...
n *= 2;
X.resize(n);