# Numerical Solutions of 2D and 3D ODE systems and phase plots

I most often use the built-in solvers in MATLAB (or the freeware alternative FreeMat) for numerical solutions of ODE systems.  They are easy to call and MATLAB graphical output is very flexible and of high quality.  That said, if I’m already working in Maxima to study stability behavior and equilibrium solutions, it’s fun to stay in the same window for  a numerical solution.

Here’s an example:  a damped, driven hard spring oscillator, written in phase space variables

$\dot{y}=z, \dot{z}=-y^3-y-z/10 + 3/10\cos(t)$

It’s straight-forward enough to run the Maxima numerical solver rk on this equation.  Extracting the numerical values of the two dependent variables is a little bit of a drag, requiring a call to makelist.  To make a plot in phase space, here’s a one-line function that is so easy to call:

/* wxphaseplot2d takes the output of s: rk([rhs1,rhs2])
and  plots in phase space */
wxphaseplot2d(s):=wxplot2d([discrete,makelist([p[2],p[3]],p,s)],[xlabel,""],[ylabel,""]);

Here’s how it works on the above  example with a cute phase-space solution curve:

And here’s another famous phase curve—the solution of the Brusselator:

And can I just point out how analytically solving for the equilibrium solution has a very similar calling sequence — so easy that there’s no reason not to do it every time you want a numerical solution, and vice versa!

Update!  Here’s a 3D phase plotter

/* wxphasplot3d takes the output of s: rk([rhs1,rhs2,rhs3]) and plots in phase space */

wxphaseplot3d(s):=wxdraw3d(point_type=0,point_size=.3,points_joined=true,points(makelist([p[2],p[3],p[4]],p,s)));

A classical example of 3D phase space trajectory curves, the Lorenz strange attractor system: