Excel macro sub srk-3()

       Sub srk-3()

       tend = Range("B5")

       step = Range("B6")

       I1 = Range("B7")

       I2 = Range("B8")

       zeta = Range("B9")

       k = Range("B10")

       T1 = Range("B11")

       T2 = Range("B12")

       c = 2# * zeta * Sqr(k * I1 * I2 / (I1 + I2))

       [F5].Select

       dt = tend / step

       ActiveCell.Offset(2, 0) = dt

       [O1].Select

       t = 0#

       x1 = 0#

       x2 = 0#

       v1 = 0#

       v2 = 0#

       For i = 1 To step

       Call rk2(dt, t, x1, x2, v1, v2, x1new, x2new, v1new, v2new, I1, I2, k, c, T1, T2)

       ActiveCell.Offset(i, 0) = t

       ActiveCell.Offset(i, 1) = x1 - x2

       ActiveCell.Offset(i, 2) = v1 - v2

       t = t + dt

       x1 = x1new

       x2 = x2new

       v1 = v1new

       v2 = v2new

       Next i

       End Sub

      Sub rk2(dt, t, x1, x2, v1, v2, x1new, x2new, v1new, v2new, I1, I2, k, c, T1, T2)

      k1 = dt * f(t, x1, v1)

      l1 = dt * g(t, x1, v1, x2, v2, I1, I2, k, c, T1)

      m1 = dt * m(t, x1, v1, x2, v2)

      n1 = dt * n(t, x1, v1, x2, v2, I1, I2, k, c, T2)

      k2 = dt * f(t + dt / 2, x1 + k1 / 2, v1 + l1 / 2)

      l2 = dt * g(t + dt / 2, x1 + k1 / 2, v1 + l1 / 2, x2 + m1 / 2, v2 + n1 / 2, I1, I2, k, c, T1)

      m2 = dt * m(t + dt / 2, x1 + k1 / 2, v1 + l1 / 2, x2 + m1 / 2, v2 + n1 / 2)

      n2 = dt * n(t + dt / 2, x1 + k1 / 2, v1 + l1 / 2, x2 + m1 / 2, v2 + n1 / 2, I1, I2, k, c, T2)

      k3 = dt * f(t + dt / 2, x1 + k2 / 2, v1 + l2 / 2)

      l3 = dt * g(t + dt / 2, x1 + k2 / 2, v1 + l2 / 2, x2 + m2 / 2, v2 + n2 / 2, I1, I2, k, c, T1)

      m3 = dt * m(t + dt / 2, x1 + k2 / 2, v1 + l2 / 2, x2 + m2 / 2, v2 + n2 / 2)

      n3 = dt * n(t + dt / 2, x1 + k2 / 2, v1 + l2 / 2, x2 + m2 / 2, v2 + n2 / 2, I1, I2, k, c, T2)

      k4 = dt * f(t + dt, x1 + k3, v1 + l3)

      l4 = dt * g(t + dt, x1 + k3, v1 + l3, x2 + m3, v2 + n3, I1, I2, k, c, T1)

      m4 = dt * m(t + dt, x1 + k3, v1 + l3, x2 + m3, v2 + n3)

      n4 = dt * n(t + dt, x1 + k3, v1 + l3, x2 + m3, v2 + n3, I1, I2, k, c, T2)

      x1new = x1 + (k1 + 2 * (k2 + k3) + k4) / 6

      v1new = v1 + (l1 + 2 * (l2 + l3) + l4) / 6

      x2new = x2 + (m1 + 2 * (m2 + m3) + m4) / 6

      v2new = v2 + (n1 + 2 * (n2 + n3) + n4) / 6

      End Sub

      Function f(t, x1, v1)

      f = v1

      End Function

      Function g(t, x1, v1, x2, v2, I1, I2, k, c, T1)

      g = k / I1 * (x2 - x1) + c / I1 * (v2 - v1) + T1 / I1

      End Function

      Function m(t, x1, v1, x2, v2)

      m = v2

      End Function

      Function n(t, x1, v1, x2, v2, I1, I2, k, c, T2)

      n = k / I2 * (x1 - x2) + c / I2 * (v1 - v2) - T2 / I2

      End Function