Write Maxima script for symbolic processing

ham.mac

@@ -0,0 +1,15 @@
+/*
+ * Copyright (c) Camden Dixie O'Brien
+ * SPDX-License-Identifier: AGPL-3.0-only
+ */
+
+ham(L, qs) := block([ps, qdots, H],
+	ps: maplist(lambda([q], concat(p_, q)), qs),
+	qdots: maplist(lambda([q], 'diff(q, t)), qs),
+	H: subst(solve(maplist(lambda([p,qdot], p = diff(L, qdot)), ps, qdots),
+	               qdots),
+	         apply("+", maplist("*", ps, qdots)) - L),
+	maplist(trigsimp, flatten(maplist(lambda([q], [
+		'diff(q,t) = diff(H, concat(p_, q)),
+		'diff(concat(p_, q),t) = -diff(H, q)
+	]), qs))));