Let's see if Sage can find a symbolic solution:
var('t')
x = function('x', t)
y = function('y', t)
u = 1.0
de1 = diff(x,t) - y == 0
de2 = diff(y,t) + x - u * y * (1 - x^2) == 0
Van_der_Pol = [de1, de2]
desolve_system(Van_der_Pol, [x, y], ics=[0, 2, 0])If you run this code, Sage will return an error. It is unable to find a symbolic solution. Now, let's try to solve it numerically:
var('t, x, y')
u = 1.0
Van_der_Pol = [y, -x + u * y * (1 - x^2)]
sol = desolve_system_rk4(Van_der_Pol, [x, y], ivar=t,
ics=[0, 2, 0], end_points=[0, 20])
t = [i for i, j, k in sol]
x_sol = [j for i, j, k in sol]
y_sol = [k for i, j, k in sol]
# Plot results
import matplotlib.pyplot as plt
plt.figure(figsize=(4, 3))
plt.plot(t, x_sol)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.subplots_adjust(bottom=0.15)
plt.savefig('Van_der_Pol_rk4.png')
plt.close()
# Limit cycle in the phase plane
plt.figure(figsize=(4, 4))
plt.plot(x_sol, y_sol)
plt.axis('scaled')
plt.xlabel('x...