Expressions for partial derivatives using AD tools

I have two bivariate functions (more generally, n n-variable functions) f_1(x_1,x_2) and f_2(x_1,x_2). I need to solve the following system of equations: \frac{\partial f_1}{\partial x_1} = 0 and \frac{\partial f_2}{\partial x_2} = 0 (the context here is a Cournot-type duopoly model). The functions are too complex to differentiate by hand, so I would like to use AD tools to solve this. However, I cannot come up with a good way to grab the partial derivatives and use them as expressions inside tools such as NLsolve or JuMP. What would be the correct way to do this?

My current approach is to calculate the Jacobian matrix, and point to the corresponding entries when solving the system. This gives me the correct solution, but feels more like a clumsy workaround than an actual solution, in particular for larger problems.

using ForwardDiff
using NLsolve

f1(x1,x2) = x1*(50-2*(x1+x2))-10-2x1
f2(x1,x2) = x2*(50-2*(x1+x2))-10-2x2

function jac(f1,f2,a)
       F(x1,x2) = [f1(x1,x2), f2(x1,x2)]
       return ForwardDiff.jacobian(x -> F(x[1],x[2]),a)

function f!(F,x)
    F[1] = jac(f1,f2,x)[1,1] # df1/dx1 == 0
    F[2] = jac(f1,f2,x)[2,2] # df2/dx2 == 0

x0 = [5.,5.] # Initial guess

res = nlsolve(f!, x0)
println( # Correct solution: [8.0,8.0]

The example above is from:

Use a closure. Eg \partial f_1 / \partial x_1 at x_1, x_2 is

ForwardDiff.derivative(x1 -> f1(x1, x2), x1)

Ah, nice and simple. Thank you. This does the trick.