How to find the regions for which the derivative of my function changes sign?

I have an energy function of the angles of the electrons’ spins. th1 is vector with (2l+2) elements and each element represents the angle of an individual electron spin. I need to eventually find the angles for which my energy is minimum. (I can use NMinimize but I want to make sure that my answer is the Global minimum, so I want to figure out the derivative first and see how many minimum my function has).

So I am taking the derivative of the function, and using Reduce to find different regions, but it takes forever ant does not give me the answer. Any idea how I can find the regions for which my derivative changes sign so I can figure out where the minimums are?

  [ScriptL]0 = 5
  [Gamma] = 
  Table[{Riffle[Range[0, -[ScriptL]0, -1], Range[[ScriptL]0]][[i]], 
  1}, {i, 1, 2 [ScriptL]0 + 1}];
  th1 = Table[Subscript[t, n] , {n, 1, 2 [ScriptL]0 + 2}]
  deriv = Table[1, {n, 1, 2 [ScriptL]0 + 2}]

  factorFxn[[ScriptL]_, m1_, m2_, p1_, p2_] := 

  If[[Gamma][[p1, 1]] - [Gamma][[m1, 1]] == [Gamma][[m2, 
 1]] - [Gamma][[p2, 1]], 
 Sum[(2 [ScriptL] + 1)^2 Sum[
 If[[Gamma][[p1, 1]] - [Gamma][[m1, 1]] == 
    mval && [Gamma][[m2, 1]] - [Gamma][[p2, 1]] == 
    mval, (-1)^([Gamma][[m1, 1]] + [Gamma][[m2, 1]] + 
      mval) ThreeJSymbol[{[ScriptL], -[Gamma][[m1, 
       1]]}, {[ScriptL], [Gamma][[p1, 
      1]]}, {[ScriptL]temp, -mval}] ThreeJSymbol[{[ScriptL]temp,
      mval}, {[ScriptL], -[Gamma][[m2, 
       1]]}, {[ScriptL], [Gamma][[p2, 
      1]]}] ThreeJSymbol[{[ScriptL], 0}, {[ScriptL], 
     0}, {[ScriptL]temp, 0}]^2, 
   0], {mval, -[ScriptL]temp, [ScriptL]temp}], {[ScriptL]temp, 
   0, 2 [ScriptL]}], 0]



 energy[th1_] :=(*(2 [ScriptL]0 +1)^2*) Sum[
 (* Find out which states we're calculating the matrix element of *)    

 (Cos[th1[[p2]]] Cos[th1[[p1]]] + 
  Cos[th1[[p2]]] Sin[th1[[p1 + 1]]] + 
  Cos[th1[[p1]]] Sin[th1[[p2 + 1]]] + 
  Sin[th1[[p2 + 1]]] Sin[th1[[p1 + 1]]] + 
  If[p1 == p2, Cos[th1[[p1]]] Sin[th1[[p1 + 1]]], 
   0]) factorFxn[[ScriptL]0, m1, m2, p1, p2]

 , {p1, 1, 2 [ScriptL]0 + 1}, {m1, 1, 2 [ScriptL]0 + 1}, {p2, 1, 
  p1}, {m2, 1, m1}];
  
 derivative = Map[D[energy[th1], #] &, th1[[1 ;; 2 [ScriptL]0 + 2]]]

  th1[[1]] = 0.000001;
  th1[[2 [ScriptL]0 + 2]] = [Pi]/2;
  Reduce[derivative == 0 && 
  0 < th1[[2 ;; 2 [ScriptL]0 + 1]] <= [Pi]/2, 
  th1[[2 ;; 2 [ScriptL]0 + 1]]]