C++ vector subscript out of range : Monotone Interpolation

Question

#include <iostream>
#include <vector>
#include <algorithm>

// Function declarations
std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq);

int main() {
    // Example usage
    std::vector<double> x = { 0,0.0208,0.0416,0.2,0.2208,0.2416,0.2416,0.2416,1.2416,1.2416,1.2624,1.2832,1.6416,1.6624,1.6832,1.6832,1.6832,2.6832,2.6832,2.704,2.7248,3.2832,3.304,3.3248,3.3248,3.3248,3.8248 };
    std::vector<double> y = { 0.404674,34.3471,31.8996,31.8996,6.74855,0.404674,0.404674,0.404674,0.404674,0.404674,30.3726,31.8996,31.8996,7.07452,0.404674,0.404674,0.404674,0.404674,0.404674,29.053,30.7901,30.7901,7.67308,0.404674,0.404674,0.404674,0.404674 };

    // Calculate xq in C++
    double dt = 0.0004;
    std::vector<double> xq;

    // Generate 10,000 equidistant samples
    double min_x = *std::min_element(x.begin(), x.end());
    double max_x = *std::max_element(x.begin(), x.end());
    for (double val = min_x; val <= max_x; val += dt)
    {
        xq.push_back(val);
    }

    std::cout << "Size of xq: " << xq.size() << std::endl;
    // Perform monotone cubic interpolation
    std::vector<double> result = wiki_mon_spline(x, y, xq);

    // Print the result
    for (size_t i = 0; i < xq.size(); ++i) {
        std::cout << "Interpolation at xq[" << i << "]: " << xq[i] << "     " << result[i] << std::endl;
    }

    return 0;
}

std::vector<double> wiki_mon_spline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& xq) {
    int n = x.size();

    // Calculate slope at grid points
    std::vector<double> h(n - 1);
    std::vector<double> delta(n - 1);

    for (int i = 0; i < n - 1; ++i) {
        h[i] = x[i + 1] - x[i];
        delta[i] = (y[i + 1] - y[i]) / h[i];
    }

    // Calculate slopes at interior points
    std::vector<double> m(n, 0.0);
    for (int i = 0; i < n - 1; ++i) {
        if (delta[i] * delta[i + 1] > 0) {
            double w1 = 2 * h[i + 1] + h[i];
            double w2 = h[i + 1] + 2 * h[i];
            m[i + 1] = (w1 + w2) / (w1 / delta[i] + w2 / delta[i + 1]);
        }
        else {
            m[i + 1] = 0.0; // Initialize to 0 if delta[i] and delta[i + 1] have opposite signs
        }
    }

    // Slope at end points
    m[0] = ((2 * h[0] + h[1]) * delta[0] - h[0] * delta[1]) / (h[0] + h[1]);

    if (m[0] * delta[0] <= 0) {
        m[0] = 0.0;
    }
    else if ((delta[0] * delta[1] <= 0) && (std::abs(m[0]) > std::abs(3 * delta[0]))) {
        m[0] = 3 * delta[0];
    }

    // Perform interpolation at query points xq
    std::vector<double> yq(xq.size(), 0.0);

    for (size_t i = 0; i < xq.size(); ++i) {
        // Find the interval containing xq[i]
        auto it = std::lower_bound(x.begin(), x.end(), xq[i]);
        if (it != x.end() && it != x.begin()) {
            int j = static_cast<int>(it - x.begin()) - 1;
            if (j < n - 1 && j >= 0) {
                // PCHIP interpolation formula
                double dx = xq[i] - x[j];
                double t = dx / h[j];
                double h00 = 2 * t * t * t - 3 * t * t + 1;
                double h10 = t * t * t - 2 * t * t + t;
                double h01 = -2 * t * t * t + 3 * t * t;
                double h11 = t * t * t - t * t;

                // Ensure monotonicity by considering signs
                double sign_delta = (delta[j] >= 0.0) ? 1.0 : -1.0;
                double sign_delta_next = (delta[j + 1] >= 0.0) ? 1.0 : -1.0;

                yq[i] = h00 * y[j] + sign_delta * h10 * h[j] + h01 * y[j + 1] + sign_delta_next * h11 * m[j];
            }
            else if (j < 0) {
                // Handle the case where xq[i] is beyond the first element of x
                yq[i] = y[0];
            }
            else {
                // Handle the case where xq[i] is beyond the last element of x
                yq[i] = y[n - 1];
            }
        }
    }

    return yq;
}

In the code above I'm performing Monotone Cubic Interpolation for 10000 samples (time 0.0004) but I'm facing error vector subscript out range however I change the limit or loop. Kindly help to resolve this? What the x input the vector range should not fall out of range that's what I expect but with a sample input x= 0 to 5 y=0 to 4 it was working fine but with original input its throwing error "Vector subscript out of range"