## Two little list utilities for MATLAB-like array indexing

I grew up with MATLAB, where extracting a subset of a vector V was easy as feeding a vector of indices into a vector.  For example entries i thru j could be had with V(i:j) and an more complicated index scheme could be accomplished with a vector of indices i_index and then V(i_index).

In Maxima, first(), last(), firstn(), rest(), and most generally makelist() allows for all that and more but with a little more cumbersome calling protocols.  Here are one-liners that achieve something like the two canonical MATLAB examples above: ## A maxima function to replicate MATLAB diff() and an efficiency comparison

I needed a function to compute a list containing the difference between adjacent entries in another list, as diff() does for an array in MATLAB.

I first wrote  something using makelist() and found it very slow for large lists.  I then rewrote the function using rest() and found it much faster: ## 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:

### wxtimeplot()

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. ### wxphaseplot2d,wxphaseplot3d

I wrote about these in another post ### sol_points()

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)

```(tw:(cosh((Vi-V3)/(2*V4))^(-1)),
m:(1/2)*(1+tanh((Vi-V1)/V2)),
w:(1/2)*(1+tanh((Vi-V3)/V4)),
aexe:aexe1*(1+tanh(Vi-V5)/V6),
ainh:ainh1*(1+tanh((Zi-V7)/V6))
)\$

(rhs1:-gCa*m*(Vi-1)-gK*Wi*(Vi-ViK)-gL*(Vi-VL)+I-ainh*Zi,
rhs2:(phi*(w-Wi)/tw),
rhs3:b*(c*I+aexe*Vi)
)\$

(V1:-.01,V2:0.15,V3:0,V4:.3,V5:0,V6:.6,V7:0,
gCa:1.1,gK:2,gL:0.5,VL:-.5,ViK:-.783,phi:.7,tw:1,
b:0.0809,c:0.22,I:0.316,aexe1:0.6899,ainh1:0.695)\$

sol:rkf45([''rhs1,''rhs2,''rhs3],[Vi,Wi,Zi],[.1,.1,.1],[t,0,400],
report=true,absolute_tolerance=1e-8)\$```
`sol_points(numsol,nth,mth):=points(map(nth,numsol),map(mth,numsol));`
```/* wxtimeplot takes the output of sol: rk45 and plots
up to 3 dependent variables
plus the scaled integration step size vs time*/
wxtimeplot(sol):=block(
[t0,t,tt,dt,big,dtbig],
t0:map(first,sol),
t:part(t0,allbut(1)),
tt:part(t0,allbut(length(t0))),
dt:t-tt,
dtbig:lmax(dt),
big:lmax(map(second,abs(sol))),
if is(equal(length(part(sol,1)),3)) then(
big: max(big,lmax(map(third,abs(sol)))),
wxdraw2d(point_type=6, key="y1",
points(makelist([p,p],p,sol)),
color=red, key="y2",
points(makelist([p,p],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel=""))
elseif is(equal(length(part(sol,1)),2)) then
wxdraw2d(point_type=6, key="y1",
points(makelist([p,p],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel="")
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",
points(makelist([p,p],p,sol)),
color=red, key="y2",
points(makelist([p,p],p,sol)),
color=green, key="y3",
points(makelist([p,p],p,sol)),
color=magenta, key="dt",point_size=.2,
points(t,big*dt/dtbig/4),xlabel="t",ylabel=""))
);```