#!/usr/bin/env python3

# integrates the van der Pol equation, making phase space and time evolution plots
# we use the Scipy ode package: "dop853" refers to an eighth-order Runge-Kutta method (basically a fancy Euler method)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
plt.rcParams.update({'font.size':14})

# DEFINE FUNCTIONS
# the model itself
def model(t,x,mu):
    xdot = np.zeros(2)
    xdot[0] = x[1]
    xdot[1] = -mu*x[1]*(x[0]**2-1)-x[0]

    return xdot

# DEFINE PARAMETERS
mu = 1

# set up start and end times of the integration
tstart = 0
tend = 100

# timestep 
dt = 0.01

# initial condition 
x0 = [1,1]

# INTEGRATE AND PLOT
t = np.linspace(tstart,tend,int((tend-tstart)/dt))
x = np.zeros((len(t),2))
x[0,:] = x0 

model_instance = ode(model).set_integrator("dop853")
model_instance.set_initial_value(x[0,:]).set_f_params(mu)

for i in range(1,len(t)):
    model_instance.integrate(model_instance.t+dt)
    x[i,:] = model_instance.y

plt.plot(t,x[:,0])
plt.xlabel('Time')
plt.ylabel('x')

plt.figure()
plt.plot(x[:,0],x[:,1])
plt.xlabel('x')
plt.ylabel('dx/dt')



plt.show()