Code: Damped Spring system

Say you have a damped spring with spring coefficient, ks =10, and damping coefficient, kd = 1 + 0.9 * sin(t). One end of the spring attaches to a fixed wall and the other end attaches to a particle. Assuming there is no other force besides the spring and damping forces applied to the particle, this is an implementation of an integration method with adaptive step size such that the particle motion is stable?

newpartpos[4] is the current ball attached to the spring
newpartpos[3] is the resting length
newpartpos[5] is the wall
kd is the damping coefficient  NOTE: the kd is unrealistic, it should be a constant. This is used to demonstrate adaptive step size.
ks is the spring coefficient

double Lr = abs(newpartpos[4][0] - newpartpos[5][0]);
double Lc = abs(newpartpos[4][0] - newpartpos[3][0]);
double kd = 1 + .9*sin(time);
double force = ks*(Lc-Lr) - velocity*kd;
double oldpartpos = newpartpos[3][0];
newpartpos[3] = newpartpos[3] + Vec3d(howFast*velocity +
                                          (0.5)*howFast*howFast*force, 0.0, 0.0);
velocity = (newpartpos[3][0]-oldpartpos)/howFast;
if (newpartpos[3][0] > newpartpos[4][0]){
velocity = -velocity;



Note that this method is the Heun method, it uses forward Euler's method as predictor and trapezoidal method as corrector. Hence, it does not need adaptive step size even though the kd is constantly changing.