Hi,

I have a series of parameterized functions containing double integrals and function limits and want to find stationary points of those functions. I have a solution using SymPy.jl, but this is only working well for k=1,2,3. For k = 4, for example, this is extremely slow (I run this for various k, vA, vB, and f(k,x)). MWE is below.

I’m wondering if I can implement this with numerical packages. I have tried with HCubature.jl but can’t get it to work, I think because the upper limit of the inner integral contains the variable of integration of the outer integral.

For differentiation I was thinking to use ForwardDiff.jl and IntervalRootFinding.jl to solve the system of first order conditions.

Any pointers how to use/connect these numerical packages would be much appreciated - or perhaps there are better solutions even. Thank you!

```
using SymPy
function func(vA,vB,wA::Sym,wB::Sym,a::Sym,b::Sym,f_a::Sym,f_b::Sym,vbar)
return integrate(integrate(f_a*f_b,(b,wB-wA+a,vbar)),(a,0,wA))*(vA-wA) +
(integrate(integrate(f_a*f_b,(b,0,wB)),(a,wA,vbar)) +
integrate(integrate(f_a*f_b,(b,0,wB-wA+a)),(a,0,wA)))*(vB-wB)
end
function findS(vA,vB,a::Sym,b::Sym,f_a::Sym,f_b::Sym,vbar,obj_func::Function)
SymPy.@vars wA wB
w = [wA, wB]
return solve(diff.(obj_func(vA,vB,w...,a,b,f_a,f_b,vbar),w), w)
end
SymPy.@vars x a b wA wB
k = 4
f(k,x) = k*x^(k-1)
println(findS(0.4,0.7,a,b,f(k,a),f(k,b),1.0,func))
```