Thank you for your response. I’ve corrected the syntax and attempted to derive the function with respect to h. This function is the general form of the one I presented previously:
f[h_,l_,d_,x_,y_]:=Sqrt[(200*Sqrt[-Cos[d]^2*Cos[h]^2*Cos[l]^2 - 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] - Sin[d]^2*Sin[l]^2 + 1]*(Cos[d]^2*Cos[h]*Cos[l]*Sin[h] + Cos[d]*Sin[d]*Sin[h]*Sin[l])*x*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2 + 10000*(Cos[d]^2*Cos[h]^2*Cos[l]^2 + 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] + Sin[d]^2*Sin[l]^2 - 1)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^3 + ((Cos[d]^2*Cos[h]^2*Cos[l]^2 + 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] + Sin[d]^2*Sin[l]^2 - 1)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^3 - (2*Cos[d]^4*Cos[h]^2*Cos[l]^2*Sin[h]^2 + 4*Cos[d]^3*Cos[h]*Cos[l]*Sin[d]*Sin[h]^2*Sin[l] + 2*Cos[d]^2*Sin[d]^2*Sin[h]^2*Sin[l]^2 - Cos[d]^2*Sin[h]^2)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]])*x^2 - ((Cos[d]^2*Cos[h]^2*Cos[l]^2 + 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] + Sin[d]^2*Sin[l]^2)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^3 - (2*Cos[d]^4*Cos[h]^2*Cos[l]^2*Sin[h]^2 + 4*Cos[d]^3*Cos[h]*Cos[l]*Sin[d]*Sin[h]^2*Sin[l] + 2*Cos[d]^2*Sin[d]^2*Sin[h]^2*Sin[l]^2 - Cos[d]^2*Sin[h]^2)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]])*y^2 - 2*(100*Sqrt[-Cos[d]^2*Cos[h]^2*Cos[l]^2 - 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] - Sin[d]^2*Sin[l]^2 + 1]*(Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l])*Sqrt[-(Cos[d]^2*Sin[h]^2 - ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2)/ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2]*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^3 + ((Cos[d]^3*Cos[h]^2*Cos[l]^2*Sin[h] + 2*Cos[d]^2*Cos[h]*Cos[l]*Sin[d]*Sin[h]*Sin[l] + Cos[d]*Sin[d]^2*Sin[h]*Sin[l]^2 - Cos[d]*Sin[h])*(-(Cos[d]^2*Sin[h]^2 - ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2)/ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2)^(3/2)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2 + (Cos[d]^5*Cos[h]^2*Cos[l]^2*Sin[h]^3 + 2*Cos[d]^4*Cos[h]*Cos[l]*Sin[d]*Sin[h]^3*Sin[l] + Cos[d]^3*Sin[d]^2*Sin[h]^3*Sin[l]^2 - Cos[d]^3*Sin[h]^3 - (3*Cos[d]^3*Cos[h]^2*Cos[l]^2*Sin[h] + 6*Cos[d]^2*Cos[h]*Cos[l]*Sin[d]*Sin[h]*Sin[l] + 3*Cos[d]*Sin[d]^2*Sin[h]*Sin[l]^2 - 2*Cos[d]*Sin[h])*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2)*Sqrt[-(Cos[d]^2*Sin[h]^2 - ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2)/ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^2])*x)*y)/((Cos[d]^2*Cos[h]^2*Cos[l]^2 + 2*Cos[d]*Cos[h]*Cos[l]*Sin[d]*Sin[l] + Sin[d]^2*Sin[l]^2 - 1)*ArcCos[Cos[d]*Cos[h]*Cos[l] + Sin[d]*Sin[l]]^3)];
solutions = Solve[{D[f[h,l,d,x,y],h]==0, -Pi < l, l < Pi, -Pi/4 < d, d < Pi/4}, h];
However, Solve and Reduce do not produce useful outputs or do not halt. Is there an automated way to detect substitution patterns, or could I use a mathematical trick to estimate the upper bound on the number of extrema for any initial values of l,d,x,y with l∈[−π,π] and d∈[−π/4,π/4]?