4th Order Adams-Moulton Method C++

Question

I am trying to solve an ODE using 4th Order Runge Kutta and the 4th Order Adams-Moulton Method. I iterate over a couple of thousand timesteps and it seems to hold fine when the values are constant(first 50 or so) However, the very first iteration when a change in input value jumps(current(I) changes from 0 to 1) it won't converge.

In order to satisfy the 4th Order Adams-Moulton Method I have split each timestep in the main function into 5 smaller timesteps so the 5th inner step will equal the next time step (0.1s) in the main function.

I will be very grateful for any advice/help!!

The code is below

#include <iostream>
#include <cmath>

using namespace std;

double state_deriv(double Vc,double I, double OCV1, double C, double Rct)
{
    double newVc = (OCV1 - Vc) / (C * Rct) - (I / C);
    return newVc;
}

double predict(double Vc,double Vc_[],double dt,double I_[],double OCV1_[],double C_[],double Rct_[])
{
    double fn[4];
    fn[0] =  state_deriv(Vc_[0], I_[0],OCV1_[0],C_[0],Rct_[0]);
    fn[1] =  state_deriv(Vc_[1], I_[1],OCV1_[1],C_[1],Rct_[1]);
    fn[2] =  state_deriv(Vc_[2], I_[2],OCV1_[2],C_[2],Rct_[2]);
    fn[3] =  state_deriv(Vc_[3], I_[3],OCV1_[3],C_[3],Rct_[3]);
    // value of next y(predicted) is returned
    double y1p = Vc_[3] + (dt/24)*(55*fn[3] - 59*fn[2] + 37 *fn[1] - 9 *fn[0]);
    return y1p;
}

double correct(double Vc,double Vc_[],double Vc1,double dt,double I_[],double OCV1_[],double C_[],double Rct_[])
{

    double e = 0.0001;
    double Vc1c = Vc1;

    double fn[4];
    fn[1] =  state_deriv(Vc_[1], I_[1],OCV1_[1],C_[1],Rct_[1]);
    fn[2] =  state_deriv(Vc_[2], I_[2],OCV1_[2],C_[2],Rct_[2]);
    fn[3] =  state_deriv(Vc_[3], I_[3],OCV1_[3],C_[3],Rct_[3]);


    do {
        Vc1 = Vc1c;
        fn[4] =  state_deriv(Vc1, I_[4],OCV1_[4],C_[4],Rct_[4]);
        Vc1c = Vc_[3] + dt/24*(9*fn[4] + 19*fn[3] - 5*fn[2] + fn[1]);
    } while (fabs(Vc1c - Vc1) > e);

    // every iteration is correcting the value
    // of state deriv using average slope
    return Vc1c;
}

double predictor_corrector(double Vc, double dt, double I[], double* OCV1, double* C, double* Rct)
{
    //double x = time[0];
    //double xn = time[1];
    double h = dt;

    //4th order RKK TO predict 4 values
    double I_[4],OCV1_[5],C_[5],Rct_[5];
    I_[0] = I[0];
    OCV1_[0] = OCV1[0];
    C_[0] = *C;
    Rct_[0] = *Rct;

    I_[4] = I[1];
    OCV1_[4] = OCV1[1];
    C_[4] = *C;
    Rct_[4] = *Rct;

    double newdt;
    for(int i = 1;i<4;i++){
        newdt = i * dt/5;
        I_[i] = I_[0] + (I[1]-I[0]) * (newdt/dt);
        OCV1_[i] = OCV1_[0] + (OCV1[1] - OCV1[0])*(newdt/dt);
        Rct_[i] = *Rct;
        C_[i] = *C;
    }

    //now do rk4 on each step
    double k1,k2,k3,k4;
    double Vc_[4];
    Vc_[0] = Vc;

    for(int i = 1;i<4;i++) {
        // 4th order runge kutta
        k1 = dt/4 * state_deriv(Vc_[i-1], *I, *OCV1, *C, *Rct);
        k2 = dt/4 * state_deriv(Vc_[i-1]+k1/2.0, I_[i], OCV1_[i], *C, *Rct);
        k3 = dt/4 * state_deriv(Vc_[i-1]+k2/2.0, I_[i], OCV1_[i], *C, *Rct);
        k4 = dt/4 * state_deriv(Vc_[i-1]+k3, I_[i], OCV1_[i], *C, *Rct);
        Vc_[i] = Vc_[i-1] + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4);
    }


    double y1p = predict(Vc,Vc_,dt/4,I_,OCV1_,C_,Rct_);
    double y1c = correct(Vc,Vc_,y1p,dt/4,I_,OCV1_,C_,Rct_);


    Vc = y1c;

    cout << Vc << endl;
    return Vc;
}

int main()
{
    double R1 = 0.0315973;
    double C = 0.00100284;
    double Vc[5];
    double OCV1[5];
    double I[5];

    Vc[0] = 4.15;
    double dt = 0.1;

    for(int i=0;i<3;i++){
        OCV1[i] = 4.15;
        I[i] = 0;
    }
    I[3]=I[4]=1;
    OCV1[3] = 4.13;
    OCV1[4] = 4.10;

    double I_[2];
    for(size_t i=1;i<4;i++){
        Vc[i] = predictor_corrector(Vc[i-1],dt,&I[i-1],&OCV1[i-1],&C,&R1);
    }
}