Plotting and solving three related ODEs in python

I have a problem that might be a little more mathematical but the crux is that I want ot solve it in python and plot it. I have three ODEs which are related to one antoher in the following way:

x''(t)=b*x'(t)+c*y'(t)+d*z'(t)+e*z(t)+f*y(t)+g*x(t)
y''(t)=q*x'(t)+h*y'(t)+i*z'(t)+p*z(t)+l*y(t)+m*x(t)
z''(t)=a*x'(t)+w*y'(t)+v*z'(t)+u*z(t)+o*y(t)+n*x(t)

How would I solve them in order to plot them in a 3d graph, through their acceleration. I know, that I have to solve them, the difficulty lies in the problem that they are second order ODEs and interlinked through their dependency onone antoher.

For some additional Info, here is the code (it doesnt really work that good feel free to try it in a different way):

from sympy import symbols, Function, Eq
import numpy as np
from scipy.integrate import solve_ivp

t = symbols('t')
x = Function('x')(t)
y = Function('y')(t)
z = Function('z')(t)

b, c, d, e, f, g = symbols('b c d e f g')
q, h, i, p, l, m = symbols('q h i p l m')
a, w, v, u, o, n = symbols('a w v u o n')

eqx = b*x.diff(t) + c*y.diff(t) + d*z.diff(t) + e*z + f*y + g*x
eqy = q*x.diff(t) + h*y.diff(t) + i*z.diff(t) + p*z + l*y + m*x
eqz = a*x.diff(t) + w*y.diff(t) + v*z.diff(t) + u*z + o*y + n*x

#First I replaced the derivatives to turn it into a first order ODE
eqx = eqx.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqy = eqy.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqz = eqz.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})

#to solve them nummerically I started to lambdify them:
freex = eqx.free_symbols
freey = eqy.free_symbols
freez = eqz.free_symbols

Xl = sp.lambdify(freex , eqx , 'numpy')
Yl = sp.lambdify(freey , eqy , 'numpy')
Zl = sp.lambdify(freez , eqz , 'numpy')

#from here on I tried to solve them, but had trouble with the dependencies and the arguments
#(so here is only the line for x)
Xo = [0]
ExampleValues = np.array([0, 0, 0, 1, 9, 0, 1, 2, 4])
space= np.linspace(0, 10, 50)
Solx = solve_ivp(XSL, (0, 10), Xo, t_eval=sace, args=ExampleValues)

Thanks for your answer in advance!

Solution

Since XSL is not defined, I can't be sure what the problem was, but here is a working version of your code (with commentary to describe the changes).

from sympy import symbols, Function, Eq
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import sympy as sp  # can't forget the imports

t = symbols('t')
x = Function('x')(t)
y = Function('y')(t)
z = Function('z')(t)
# these symbols needed to be defined before use
dx, dy, dz = symbols('dx, dy, dz')

# Unchanged
b, c, d, e, f, g = symbols('b c d e f g')
q, h, i, p, l, m = symbols('q h i p l m')
a, w, v, u, o, n = symbols('a w v u o n')

# Unchanged
eqx = b*x.diff(t) + c*y.diff(t) + d*z.diff(t) + e*z + f*y + g*x
eqy = q*x.diff(t) + h*y.diff(t) + i*z.diff(t) + p*z + l*y + m*x
eqz = a*x.diff(t) + w*y.diff(t) + v*z.diff(t) + u*z + o*y + n*x

# Unchanged
eqx = eqx.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqy = eqy.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})
eqz = eqz.subs({ sp.Derivative(y,t):dy, sp.Derivative(z,t):dz, sp.Derivative(x,t):dx,})

# Here is where the main changes start
# To keep things simple, pass all variables to all lambda functions
# The first six are the state variables (the original variables and 
# their derivatives, since this is second order).
# The remainder are parameters.
free = [x, y, z, dx, dy, dz, 
        b, c, d, e, f, g, q, h, i, p, l, m, a, w, v, u, o, n]

# You might have forgotten these equations? They define the (trivial)
# relationship between
# the derivatives of the first three state variables (`x`, `y`, `z`)
# and the last three state variables (`dx`, `dy`, `dz`).
dXl = sp.lambdify(free, dx , 'numpy')
dYl = sp.lambdify(free, dy , 'numpy')
dZl = sp.lambdify(free, dz , 'numpy')

# I've just prepended `dd` to the names of these to indicate that 
# they are expressions for the second time derivatives.
ddXl = sp.lambdify(free, eqx , 'numpy')
ddYl = sp.lambdify(free, eqy , 'numpy')
ddZl = sp.lambdify(free, eqz , 'numpy')

# The "right hand side" function
# Accepts the state variables passed by `solve_ivp` as a vector
# Accepts the parameters individidually
# Evaluates and returns the derivatives of the state variables
# as defined above
def df(t, state, *params):
    return [dXl(*state, *params), dYl(*state, *params), dZl(*state, *params),
            ddXl(*state, *params), ddYl(*state, *params), ddZl(*state, *params)]

# Here I've set values to zero so that we have three identical,
# *independent* spring-mass-damper systems, just so I know the
# behavior to expect. You can set them as you wish to couple the 
# systems.
c = d = q = i = a = w = 0
f = e = m = p = n = o = 0

# Negative to get oscillatory solution
# Alternatively, include the negative signs somewhere above
b = h = v = -0.1  # damping
g = l = u = -1  # mass

# Be sure to keep the order consistent with `free` above
x = y = z = 1
dx = dy = dz = 0
Xo = [x, y, z, dx, dy, dz]
params = b, c, d, e, f, g, q, h, i, p, l, m, a, w, v, u, o, n

t_eval = np.linspace(0, 10, 50)
sol = solve_ivp(df, (0, 10), Xo, t_eval=t_eval, args=params)

t = sol.t
x, y, z, dx, dy, dz = sol.y
plt.plot(t, x)

(Note that that the use of Function and Derivative were not really needed here; this would work just as well defining the second derivatives in terms of regular symbols x, y, z, dx, dy, and dz. I've left that part alone, though. Perhaps you needed to start with Functions if you are deriving the real equations with SymPy.)