Using Sympy to solve System of ODE

I would like to model a double pendulum with Lagrangian mechanics. The angles phi1 and phi2 are the variables in my model. Application of Lagrange leads to a system of ordinary equations: Two equations with the variables phi1'', phi1', phi1 and phi2'', phi2', phi2.

So far I've added placeholder variables for all derivatives, substituted the current values and solved for the second derivatives (the only unknowns remaining). Then I used explicit euler. Here's my result: https://www.youtube.com/watch?v=Qh-faLBYhBA&feature=youtu.be

Instead of explicit euler, I would like to use a higher order method.

Sympy has the dsolve function. Looks very nice, but I can't figure out how to use it. You pass a list of equations and a list of functions, but this leads to a cryptic error message:http://nbviewer.jupyter.org/gist/lhk/e58bc62256a9d0dbad188c922cb2a64a

What is going wrong ?