I do know that we can differentiate equations on our own just using R.

Load expression

(my_expression <- expression(x ^ 2 * y + 3 * x ^ 2 + 4 * y - 5))
## expression(x^2 * y + 3 * x^2 + 4 * y - 5)

Differentiate

D(my_expression, 'x')
## 2 * x * y + 3 * (2 * x)

Differentiate twice, different variables:

D(D(my_expression, 'x'), 'y')
## 2 * x

On intital inspection of the cubature package, it only seems that this libary integrates over multiple dimensions, not differentiate (http://homepages.math.uic.edu/~jyang06/stat401/handouts/handout8.pdf), so i’m not too sure how to use it for this specific problem. However, this is a really interesting package. Thanks for pointing me to it and i’ll look more into it when I complete the final project submittal.