如何在Python中绘制耦合的二阶非线性常微分方程的二阶导数？

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def equations(t, y):
    x1, x2, v1, v2 = y
    dx1dt = v1
    dx2dt = v2
    dv1dt = -x1**3 - x2
    dv2dt = -x1 - x2**3
    return [dx1dt, dx2dt, dv1dt, dv2dt]

y0 = [1, 2, 0, 0]  # 初始条件
t_span = [0, 10]  # 时间范围

sol = solve_ivp(equations, t_span, y0)

plt.plot(sol.t, sol.y[0], label='x1')
plt.plot(sol.t, sol.y[1], label='x2')
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.show()