I am dealing with the following piece of code, to study a problem in quantum mechanics:

`L[n_] := KirchhoffMatrix[CompleteGraph[n]]; c[n_] := 1; w[n_, p_] := Table[KroneckerDelta[k, p], {k, 1, n}]; P[n_, p_] := KroneckerProduct[w[n, p], w[n, p]]; s[n_] := 1/Sqrt[n]*Table[1, {k, 1, n}]; Ps[n_] := KroneckerProduct[s[n], s[n]]; H[n_, \[Lambda]_] := \[Lambda] L[n] - P[n, c[n]]; U[n_, t_, \[Lambda]_] := MatrixExp[-I*t*H[n, \[Lambda]]]; \[Psi][n_, t_, \[Lambda]_] := U[n, t, \[Lambda]].s[n]; prs[n_, t_, \[Lambda]_] := Abs[w[n, c[n]].\[Psi][n, t, \[Lambda]]]^2; Prob[n_] := NMaximize[prs[n, t, \[Lambda]], {t,\[Lambda]}][[1]] `

The NMaximize function takes quite a while on my machine to compute $ \text{Prob}(n)$ , so I would be interested in any suggestion to increase the efficiency of the code – taking into account that the input graph could be a different one. Probably the hardest part to compute is taking the exponential matrix of $ H$ , but I’m not aware of any way to optimize it.