A general solution can be obtained simply by choosing different sets of signs in a quaternion norm identity, analogous to how the general solution of the sum of two squares equation of the same type can be obtained by considering norms of Gaussian integers.

Many thanks for your answer Carlo (which I've marked as the accepted one), and your kind offer to email me this paper. That would be great, & my email address is [email protected]

Well both, i.e. discrete space and time, if that makes sense for a dynamic system. Obviously the dynamic laws wouldn't involve derivatives, but then that's what I'm asking about.