Can you specify a bit more precisely exactly what you need? Here’s a starter example with TaylorSeries.jl (which is a registered package, so just Pkg.add("TaylorSeries") to install):
julia> using TaylorSeries
julia> f(x, y) = cos(x*y) + sin(x - y^2);
julia> xx, yy = set_variables("x y", order=15);
julia> f(xx, yy)
1.0 + 1.0 x² - 1.0 y² - 0.5 x² y² - 0.16666666666666666 x⁶ + 0.5 x⁴ y² - 0.5 x² y⁴ + 0.16666666666666666 y⁶ + 0.041666666666666664 x⁴ y⁴ + 0.008333333333333333 x¹⁰ - 0.041666666666666664 x⁸ y² + 0.08333333333333333 x⁶ y⁴ - 0.08333333333333333 x⁴ y⁶ + 0.041666666666666664 x² y⁸ - 0.008333333333333333 y¹⁰ - 0.0013888888888888887 x⁶ y⁶ - 0.00019841269841269839 x¹⁴ + 0.0013888888888888885 x¹² y² - 0.004166666666666666 x¹⁰ y⁴ + 0.006944444444444443 x⁸ y⁶ - 0.006944444444444443 x⁶ y⁸ + 0.004166666666666666 x⁴ y¹⁰ - 0.0013888888888888885 x² y¹² + 0.00019841269841269839 y¹⁴ + 𝒪(‖x‖¹⁶)
There are various derivative methods available to extract derivatives, or you can just access the coefficients of this multivariate polynomial.
If you need to expand e.g. around (-3, 4), you can do
f(3+xx, -4+yy)
Here, xx and yy now represent small perturbations around -3 and 4, respectively.