Thank you for your reply. I made a (in any case irrelevant) mistake in the definition of omsq[k]. I report the correct code:
m = 10;
ll = 1;
mro = 1;
mpi = 1;
g = 5.5;
vl := Table[{i, j, l}, {i, -m, m}, {j, -m, m}, {l, -m, m}];
vll := Flatten[vl, 2];
s[v_, k_] := 1/(v^2 - (k*ll/(2*Pi))^2);
ss[k_] := Total[Map[s[Norm[#], k] &, vll]]
phi[k_] = ArcTan[-2*(Pi)^2*k*ll/(2*Pi*ss[k])];
omsq[k_] := ((mro)^2 + k^2)*4;
gamma[k_] := g^2*k^3/(6*Pi*omsq[k]);
h1[k_] :=
g^2*k^3*Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/(3*(Pi)^2*
Sqrt[omsq[k]]);
h2[k_] :=
g^2*k^2*(1 + (1 + 2*(mpi)^2/omsq[k])*Sqrt[omsq[k]]*
Log[(Sqrt[omsq[k]] + 2*k)/(2*mpi)]/k)/(6*(Pi)^2*
Sqrt[omsq[k]]);
a0 := h1[mro] - mro*h2[mro]/2 + g^2*(mpi)^2/(6*(Pi)^2)
a[k_] :=
h1[mro] + (omsq[k] - (mro)^2)*h2[mro]/(2*mro) - h1[k] +
I*Sqrt[omsq[k]]*gamma[k];
fpi[k_] := ((mro)^2 - a0)/((mro)^2 - omsq[k] - a[k])
delta[k_] :=
ArcCot[((mro)^2 - omsq[k] -
h1[mro] - (omsq[k] - (mro)^2)*h2[mro]/(2*mro) + h1[k])/(omsq[k]*
gamma[k])];
However, I'm looking for an instruction that automatically finds a solution without the need to set a starting point, so NSolve seems to be my choice. And in fact, I already tried it before opening this question, but the point is that if I write
NSolve[delta[k] + phi[k] == 0 && k >= 0, k, Reals]
the process seems to be too long. I waited 5 or 10 minutes whereupon I aborted the process. Looking at the plot of the two functions, I see that all the solutions are in the range k<=120. If, instead I write
NSolve[delta[k] + phi[k] == 0 && 0 <= k < 120, k, Reals]
similar to how you suggested, it finds in 1 minute all these solutions:
{{k -> 0.}, {k -> 60.2778}, {k -> 60.7}, {k -> 61.2595}, {k ->
61.6217}, {k -> 61.9575}, {k -> 62.3592}, {k -> 62.681}, {k ->
62.8904}, {k -> 63.3799}, {k -> 63.7477}, {k -> 64.1951}, {k ->
64.5501}, {k -> 64.8526}, {k -> 65.1896}, {k -> 65.4012}, {k ->
65.7166}, {k -> 66.4009}, {k -> 66.93}, {k -> 67.2742}, {k ->
67.5128}, {k -> 67.8014}, {k -> 68.1651}, {k -> 68.6514}, {k ->
68.9567}, {k -> 69.2618}, {k -> 69.6373}, {k -> 69.9291}, {k ->
70.4181}, {k -> 70.9846}, {k -> 71.1348}, {k -> 71.5875}, {k ->
71.7683}, {k -> 72.0938}, {k -> 72.349}, {k -> 72.5804}, {k ->
73.1449}, {k -> 73.418}, {k -> 73.7302}, {k -> 73.9907}, {k ->
74.2564}, {k -> 74.5236}, {k -> 74.7941}, {k -> 75.2711}, {k ->
75.4883}, {k -> 75.7975}, {k -> 76.1651}, {k -> 76.4205}, {k ->
76.8743}, {k -> 77.3373}, {k -> 77.5906}, {k -> 77.9162}, {k ->
78.1684}, {k -> 78.8594}, {k -> 79.4396}, {k -> 79.9125}, {k ->
80.1565}, {k -> 80.3452}, {k -> 80.5801}, {k -> 80.8826}, {k ->
81.3629}, {k -> 81.8248}, {k -> 82.1206}, {k -> 82.3513}, {k ->
82.5568}, {k -> 82.8177}, {k -> 83.5153}, {k -> 83.7607}, {k ->
83.9889}, {k -> 84.1856}, {k -> 84.4299}, {k -> 84.7101}, {k ->
85.2929}, {k -> 85.8794}, {k -> 86.2753}, {k -> 86.5458}, {k ->
87.0464}, {k -> 87.4104}, {k -> 88.0525}, {k -> 88.3412}, {k ->
88.7546}, {k -> 89.0501}, {k -> 89.6893}, {k -> 90.0861}, {k ->
90.7352}, {k -> 91.2202}, {k -> 91.6101}, {k -> 92.2823}, {k ->
92.4872}, {k -> 94.1707}, {k -> 94.4087}, {k -> 94.8186}, {k ->
95.1993}, {k -> 96.4693}, {k -> 97.9171}, {k -> 98.2287}, {k ->
99.0925}, {k -> 101.624}, {k -> 102.032}, {k -> 105.249}, {k ->
108.799}}
but it missed all the solutions for k<60 !!!
I don't know why every time I insert an interval for k in which there are very many solutions, Mathematica misses all the ones in the first part of the interval. I suppose this problem can be overcome with special instruction. Maybe I have to request a finer calculation? I would insert such a command, but at the same time, I don't need all the solutions. So, to save time, I would ask Mathematica to find only the first, say, ten solutions. In this way, I don't even need to introduce an "ad hoc" upper limit for k.
How can I do this? Thank you.
