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)
end
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
end
x0 = [5.,5.] # Initial guess
res = nlsolve(f!, x0)
println(res.zero) # Correct solution: [8.0,8.0]
```

The example above is from: http://faculty.econ.ucdavis.edu/faculty/bonanno/teaching/200C/Cournot.pdf