Como solucionar sistema de ecuaciones diferenciales complejas numéricamente

Tengo un sistema de dos ecuaciones diferenciales complejas no lineales, me gustaría saber el comportamiento de la solución, para ello opté por intentar resolverlas numéricamente, pero me encontré con el obstáculo que no sabía programar sistemas de ecuaciones complejas.
Entonces se me ocurrió separar la parte real de la imaginaria de la compleja obteniendo una sistema de 4 ecuaciones diferenciales ordinarias no lineales reales, intenté resolver este sistema por runge-kutta de orden 4 pero no he podido implementarlo de la manera adecuada y además la manipulación algebraica de la separación de las partes imaginarias y relaes del sistema original fue bastante engorroso y largo por lo cual es muy propenso a errores.
Me gustaría saber si alguién conoce alguna forma de resolver un sistema de ecuaciones ordinarias no lineales complejas de manera directa en algún lenguaje de programación o paquetería.
Por lo que he leido en pylab, con el comando "odeint" es posible, pero sólo he podido resolver ecuaciones ordinarias de primer orden.