I would recommend that you use the Scipy module for this. Your equations seem to be a simple set of linear ordinary differential equations. You'll have to define a straightforward function in python to return the 8 derivatives (that is, p'1(t)...p'8(t)). Then, as your system defines an initial value problem, you can use the solve_ivp function from scipy.integrate. You'll have to pass the derivative outputting function defined by you, the span of time over which you want a solution and the initial condition. The solve_ivp will give you the solution.
Here's an excellent tutorial from the Scipy developers on how to use solve_ivp for the Lorenz equations. Try to walk through the tutorial and replicate the results. After that, you just need to replace the Lorenz equations in that code with your own...and your problem is solved!
Joakim Subdnes wrote an excellent note on exactly this, with a nice example (starting on Page 31) of solving systems with only numpy. See https://sundnes.github.io/solving_odes_in_python/ode_book.pdf
If you're case doesn't need to be terribly general the SymPy library can be very helpful (https://docs.sympy.org/latest/modules/solvers/ode.html)