Runge-Kutta Behaving Strangely When Initial Value Is Negative

I have the following RK4 algorithm in C++ for solving first-order differential equations:

using ODE_Function = std::function<double/*dy/dx*/(double/*x*/, double/*y*/)>;

template<int length> //"length" is the number of data points
double* rk4(ODE_Function fxn, double y0, double x0, double h) {

    double* y = new double[length];
    
    y[0] = y0; //initialize output array
    double x = x0; //n is our discretized x variable

    //main loop
    for (int n = 1; n < length; n++) {//y0 is already known, so we start from the second data point
        
        double k1 = fxn(x, y[n-1]);
        double k2 = fxn(x + h/2, y[n-1] + h*k1/2);
        double k3 = fxn(x + h/2, y[n-1] + h*k2/2);
        double k4 = fxn(x + h, y[n-1] + h*k3);

        y[n] = y[n-1] + h/6*(k1 + 2*k2 + 2*k3 + k4);
        x += h;
    }

    return y;
}

However, this behaves strangely when x0 is negative...

Specifically, the following cases should theoretically output results that follow the trend of (-x^2)-2(x+1), just viewed from different "windows". But only the final one does this. The first one results in a negative exponential, and the second one results in a positive exponential for some reason.

double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -1; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}

double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -5; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}

double fxn(double x, double y) {
    return x*x + y;
}
int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = 0; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);
}

There is nothing wrong with your Runge-Kutta implementation - just your expectation of the solution.

Rearrange the equation slightly and use an integrating factor. You will find an exact solution

y = (y₀+x₀²+2x₀+2)exp(x-x₀)-x²-2x-2

Try it with your main() function. (You will need to include iostream, functional and cmath headers.)

int main() {
    const int length = 100;
    double y0 = -2;
    double x0 = -5; //LOOK HERE
    double h = 0.1;

    double* result = rk4<length>(fxn, y0, x0, h);

    for ( int i = 0; i < length; i++ ) 
    {
       double x = x0 + i * h;
       double exact = ( y0 + x0 * x0 + 2 * x0 + 2 ) * exp( x - x0 ) - x * x - 2 * x - 2;
       cout << i << '\t' << x << '\t' << result[i] << '\t' << exact << '\n';
    }
}