$$ \ddot{x}_1 = f_1(x_1,\dot{x}_1,x_2,\dot{x}_2) $$ $$ \ddot{x}_2 = f_2(x_1,\dot{x}_1,x_2,\dot{x}_2,x_3,\dot{x}_3) $$ $$ \ddot{x}_3 = f_3(x_2,\dot{x}_2,x_3,\dot{x}_3,x_4) $$ $$ \ddot{x}_4 = f_4(x_3,x_4,\dot{x}_4,x_5,\dot{x}_5) $$ $$ \ddot{x}_5 = f_5(x_4,\dot{x}_4,x_5,\dot{x}_5,x_6) $$ $$ \ddot{x}_6 = f_6(x_5,x_6) $$
I've got the above stiff system of 2nd order differential equations and I know all of the initial conditions. I solve this by converting the system into first order differential equations and then using numerical methods I integrate (I've tried both forward and backward Euler and explicit RK4). At steady-state everything is fine but one of the constants in $f_4$ and $f_5$ make an abrupt change during the simulation and I'm not sure how deal with it.
I've tried making the step size very small but the system responds wildly with all three methods. Is my solution procedure incorrect?