Implementing the Lanczos algorithm into C++ for a quantum anharmonic oscillator

Question

Firstly, I would like to mention that I am a complete beginner when it comes to coding, let alone C++, so bear with me, as I need complete guidance. My task is to implement the Lanczos algorithm for the case of a 1-D anharmonic oscillator in C++, with reference to the paper linked Analytical Lanczos method.

The paper offers a step by step guide for the implementation of the algorithm: Step by step guide here

with the initial trial function being: Psi_1 = (1 + x^2) * (exp(-x^2 - 1/4 * x^4).

The paper also contains code in MATHEMATICA for this particular case. Mathematica code

and thus, here is my attempt, which is greatly unfinished, however, I wanted to ensure I was going along the correct path with regards to the programming logic. There are still plentiful errors etc. (Also excuse the lack of fundamentals here, I am only a beginner. Thank you very much.)

int main() {

    //Grid parameters.
    const int Rmin = 1, Rmax = 31, nx = 300;//Grid length and stepsize.

    double dx = (Rmax- Rmin) / nx; //Delta x.

    double a, b;

    std::vector<double> x, psi_1;

    for (int j = 1; j < 64; ++j) { //Corresponds to each succesive Lanczos Vector.

    
        for (int i = Rmin; i < nx + 1; i++) { //Defining the Hamiltonian on the grid. 

            x[i] = (nx / 2) + i;

            psi_1[i] = (1 + pow(x[i] * dx, 2)) * exp(pow(-x[i] * dx, 2) - (1 / 4 * pow(x[i] * dx, 4 )) //Trial wavefunction.

            H[i] = ((PSI[j][i + 1] - 2 * PSI[j][i] + PSI[j][i - 1]) / pow(dx, 2)) + PSI[j][i] * 1/2 * pow(x[i] * dx, 2) + PSI[j][i] * 2 * pow(x[i] * dx, 4) + PSI[j][i] * 1/2 * pow(x[i], 6); //Hamiltonian. ****   

            //First Lanczos step.

            PSI[1][i] = psi_1[i]
        }

    //Normalisation of the wavefunction (b).

    double b[j] = 0.0; 

    for (int i = Rmin; i < nx + 1; i++) { 

        PSI[1][i] = psi_1[i];

        b[j] += abs(pow(PSI[j][i], 2));
    }

    b[j] = b[j] * dx;

    for (int i = Rmin; i < nx + 1; i++) {

        PSI[j] = PSI[j] / sqrt(b[j]);
    }
    //Expectation values (a). Main diagonal of the Hamiltonian matrix.

    double a[j] = 0.0;

    for (int i = Rmin; i < nx + 1; i++) {

        a[j] += PSI[j] * H[i] * PSI[j] * dx
    }
    
    //Recursive expression.

    PSI[j] = H[i] * PSI[j-1] - PSI[j-1] * a[j-1] - PSI[j-2] * b[j-1]

    //Lanczos Matrix.

    LanczosMatrix[R][C] = 
    for (int R = 1; R < 64; R++) {
        row[R] = 
    }
}  

I have yet to finish the code, but some experienced guidance would be greatly appreciated! (also, the code has to be cleaned up greatly, but this was an attempt to get the general idea down first.)