In divround, pint and qint must be the integer parts of pfrac and qfrac, However, as I understand, in Python, converting a float to an int always rounds it towards 0, i.e., positive floats are rounded down, but negative floats are rounded up, so int(-2.7) is -2, not -3. Instead of
pint, qint = [int(thing) for thing in (pfrac, qfrac)]
you should have
A meshless numerical library C++ Medusa is claimed to be dimension indenpendent. The implementation itself strongly depends on the C++ templates.
For example in this paper they use Medusa to sole a Poisson problem in 1, 2, 3, and 4 dimensional space, with basically the same code.
I suspect it might be easier to only generate valid expressions if you start with an infix expression "tree" of your expression and convert that to RPN format if you really want to.
Each node can either be a binary operator $O$ or a number $N$. $O$ nodes by definition cannot be leaf nodes since they need exactly two child nodes, and $N$ nodes must ...
Your code has a small issue that is distorting the amplitude. You are not initializing p, which throws off the Verlet steps right from the beginning. If you have this set, the amplitude will stay roughly between 1 and -1.
import matplotlib.pyplot as plt
import numpy as np
e = 200
dx = 0.08
xlim = 100
p = ...