Plotting Numerical Solutions of Ordinary Differential Equations in Maxima

Maxima comes with a few well-implemented numerical ODE solvers — rk() and rkf45() — and I’ve written in a previous post about my variable stepsize implementation of the 2nd order backward differentiaion formula for stiff systems  BDF2a().

Plotting the resulting numerical solutions — which is a list of lists of the form

[[indvar, depvar1, depvar2,… ], [indvar, depvar1, depvar2,… ], …   ]

with each of the sublists  representing the state of the system at a discrete time — requires a little familiarity with list functions like makelist() and map().

Here are a few little helper routines to make that quicker:


This is mainly a teaching tool for my ordinary differential equations course, taking the output of the solver and, assuming there aren’t more than three dependent variables, plotting all the dependent variables vs the independent variable on the same set of axes, with a representation of the step size.



I wrote about these in another post



This is a more general purpose function.  Given the output of the numerical ode solver, this constructs a points() subcommand ready to plug into draw2d().  This allows for the full generality of draw2d:  colors, axis labels, point points,

Especially useful for plotting results from more than one call to the solver, such as when we need to see the effect of changing a parameter value:


Below I’ve pasted the code for the plots, and also for the interesting differential equations from  Larter et al in Chaos, V.9 n.3 (1999)





/* wxtimeplot takes the output of sol: rk45 and plots 
 up to 3 dependent variables 
plus the scaled integration step size vs time*/
if is(equal(length(part(sol,1)),3)) then(
 big: max(big,lmax(map(third,abs(sol)))),
 wxdraw2d(point_type=6, key="y1",
 color=red, key="y2",
 color=magenta, key="dt",point_size=.2,
elseif is(equal(length(part(sol,1)),2)) then
 wxdraw2d(point_type=6, key="y1",
 color=magenta, key="dt",point_size=.2,
elseif is(equal(length(part(sol,1)),4)) then(
 big: max(big,lmax(map(third,abs(sol)))),
 big: max(big,lmax(map(fourth,abs(sol)))),
 wxdraw2d(point_type=6, key="y1",
 color=red, key="y2",
 color=green, key="y3",
 color=magenta, key="dt",point_size=.2,

BDF2a: An adaptive-stepsize Stiff ODE solver in Maxima

**March 2021 update** Maxima has a full-featured stiff solver called Eulix now, thanks to Helmut Jarausch.


Since 2011, Maxima has included the user-contributed numerical ODE solver rkf45() created by Panagiotis Papasotiriou.  This implementation of the fourth (and fifth) order Runge-Kutta-Fehlberg embedded method features adaptive timestep selection and a nicely optimized function evaluation to make it run pretty fast in Maxima.  As an explicit RK method, it is suitable for non-stiff equations.  To use it in Maxima, you must first load(rkf45).

For years in my differential equations class, we’ve test-driven the various numerical ODE methods in MATLAB and explored their behavior on stiff systems.  This year we’ve switched from MATLAB to Maxima, and I’ve written a stiff solver for that purpose.  It’s an adaptive stepsize  implementation of the second-order Backward Difference Formula.  The BDF.mac file can be found here.  The package contains a fixed-stepsize implementation BDF2()  and the adaptive stepsize BDF2a().

I’ve tried to make the calling sequence and options as close to rkf45() as possible.  I have not made any attempt to optimize the performance of the method.  BDF methods are implicit, meaning that a nonlinear equation (or system of equations) must be solved at each step.  For that I use the user-contributed Maxima function mnewton().

Additionally, note that I included my homemade help utility so that once BDF.mac is loaded, help(BDF2) and help(BDF2a) display help lines in the console window. For simplicity in comparing the methods, the load(rkf45) command is included in BDF.mac, as is an rkf45 help file for use with help().   Finally, I included in the package my little solution plotting utility wxtimeplot() which takes the “list of lists” output of BDF2(), BDF2a() and rkf45() and plots the dependent variables vs the independent variable along with a log-scale plot of the stepsize.

Below is an example that illustrates how the the BDF method performs relative to rkf45 for an increasingly stiff problem.

The equation is

 \dot{y} = -\lambda (y-\sin (10t)-t)+10 \cos (10t)+1 .    

For small values of  \lambda , the explicit rkf method’s adaptive stepsizes result in much more efficient solutions for a given accuracy level.  However for large \lambda the implicit BDF method is comparatively much more efficient, in terms of number of steps, and  number of function evaluations, although in its current implementation, the BDF method takes more time.

Here’s the results for \lambda=50:


But for \lambda=5000, the performance situation is reversed, with the BDF formula essentially immune from the effects of extra stiffness, while the rkf method requires an order of magnitude more steps.



E. Alberdi Celaya , J. J. Anza Aguirrezabala , and P. Chatzipantelidis.  “Implementation of an adaptive BDF2 formula and comparison with the MATLAB Ode15s”,  Procedia Computer Science,  December 2014