How to get the state of a dynamic system from lyapunovspectrum?

General Usage
1

Hello everyone, I am writing a code for calculating the Lyapunov indicators of the system:

ks = 0.0:0.0001:1.0; λs = zeros(length(ks), 6)
for (i, k) in enumerate(ks)
    set_parameter!(ds, 14, k)
    λs[i, :] .= lyapunovspectrum(ds, 300000; Δt = 0.01, show_progress = true)
end

As far as I know, the lyapunovspectrum function integrates the system each time for each new value of the parameter k. So, I want to make sure that with each new integration of the system, the state of the system at the last step of the previous integration is taken as the initial conditions. In python, I wrote this:

initial_state = np.array([-1.5, 0.0, 0.0, -2.5, 0.0, 0.0])

for j, k in enumerate(ks):

  twoHRs = [
          y(1) - a * y(0) ** 3 + b * y(0) ** 2 - y(2) + I - k1 * (y(0) - v_s) * (1 / (1 + sy.exp(-lambd * (y(3) - Theta)))) + k * (y(3) - y(0)),
          c - d * y(0) ** 2 - y(1),
          r * (s * (y(0) - x_r) - y(2)),
          y(4) - a * y(3) ** 3 + b * y(3) ** 2 - y(5) + I - k2 * (y(5) - v_s) * (1 / (1 + sy.exp(-lambd * (y(0) - Theta)))) + k * (y(0) - y(3)),
          c - d * y(3) ** 2 - y(4),
          r * (s * (y(3) - x_r) - y(5))
          ]
  
  print('k = ', k)
  n = len(twoHRs)
  ODE = jitcode_lyap(twoHRs, n_lyap=n)
  ODE.set_integrator("vode")
  ODE.set_initial_value(initial_state,0.0)
  times = range(1,300000,1)
  lyaps = []
  for time in times:
    lyaps.append(ODE.integrate(time)[1])
    # converting to NumPy array for easier handling
  lyaps = np.vstack(lyaps)

  #we take values from the last point for new initial conditions
  final_state = np.array([ODE.y_dict.get(y(0)), ODE.y_dict.get(y(1)), ODE.y_dict.get(y(2)), ODE.y_dict.get(y(3)), 
                        ODE.y_dict.get(y(4)), ODE.y_dict.get(y(5))])
  
  initial_state = final_state
  
  temp = np.zeros(6)

  for i in range(n):
    lyap = np.average(lyaps[1000:,i])
    temp[i] = lyap
    stderr = sem(lyaps[1000:,i]) # Note that this only an estimate
    print("%i. Lyapunov exponent: % .4f ± %.4f" % (i+1,lyap,stderr))
  
  lce1[j] = temp[0]
  lce2[j] = temp[1]
  lce3[j] = temp[2]
  lce4[j] = temp[3]
  lce5[j] = temp[4]
  lce6[j] = temp[5]

Tell me, please, is it possible to do the same on Julia?