(Hyper) Dual Numbers for automatic differentiation (CAS)

[mcasl]

Dear all,
I'm new to HP Prime and I am taking my first steps in programming it.
In the context of CAS, I'm working on a way of automating the automatic differentiation of a function using dual numbers.
For example, let \[f\left(x,y\right)=\left(x^2+7y\right)*\left(7x^2+5y\right)\] . To find the derivative at x=3, y=2 we apply the function at \(x=3+\epsilon_1\) and \(y=2+\epsilon_2\), giving, after simplification:
\[f\left(3+\epsilon_1,2+\epsilon_2\right)=7\epsilon_1^4+84\epsilon_1^3+54\epsilon_​1^2 \epsilon_2+486 \epsilon_1^2+324\epsilon_1\epsilon_2+1404\epsilon_1+35\epsilon_2+626\epsilon_2+1​679 \]
Now, that expression can be further simplified considering that, by definition, \(\epsilon_1^2=\epsilon_2^2=\epsilon_1\epsilon_2=0\), yielding:
\[f\left(3+\epsilon_1,2+\epsilon_2\right)=1679+1404\epsilon_1+626\epsilon_2\]
The result shows that \(f\left(3,2\right)=1679\) and that \(\frac{\partial f}{\partial x}=1404\) and \(\frac{\partial f}{\partial y}=626\).
I'm currently using the CAS function taylor2d that I downloaded from the forums in order to discard the terms with \(\epsilon_i^n |_{ n>2}\) as well as \(\epsilon_i\epsilon_j\), but that approach does not cover multivariate cases with more than two variables.
I'm interested in knowing if there is an alternative way of simplifying the result, especially if it involves letting xcas know that the rules for the \(\epsilon_i\) are the ones previously mentioned.
Thanks for your patience and attention