Lieb-Liniger model and Bethe ansatz

Solver for the Bethe ansatz equation

(*c:Interaction strength*)
c=5.0;
(*Bethe roots: Num:# particles, L:length of the ring, Int:integer set*)
BetheRoots[Num_,L_,Int_]:=Table[kset[i],{i,1,Num}]/.Module[{Iset},
Iset[i_]:=Int[[i]];
FindRoot[Table[kset[j] L==2 π Iset[j]-2 Sum[ArcTan[(kset[j]-kset[l])/c],{l,1,Num}],{j,1,Num}],
Table[{kset[j],(-(Num-1)/2+(j-1)) 2 π/L},{j,1,Num}]
]
]
(*Dispersion relation*)
Disp[Num_,L_,Int_]:=Module[{root},
root:=BetheRoots[Num,L,Int];
{Total[root],Total[root^2]}
]

Dispersion relation (Low-lying states)

Nump=16(*=Number of particles. The length of the ring = Nump*);
IntList=Module[{lis,inilis,dim},
lis=PadLeft[#,Nump] &/@Map[Sort,Flatten[Table[IntegerPartitions[i],{i,0,Nump}],1]]//Parallelize;
inilis=Table[-(Nump-1)/2+(j-1),{j,1,Nump}];
dim=Length[lis];
Table[inilis+lis[[i]],{i,1,dim}]
];
DispTab=Module[{pretab,dim,minE},
pretab=Disp[Nump,Nump,#] &/@IntList//Parallelize;
dim=Length[pretab];
minE=Min[Table[pretab[[k]][[2]],{k,1,dim}]];
Join[Table[{pretab[[k]][[1]],pretab[[k]][[2]]-minE},{k,1,dim}],Table[{-pretab[[k]][[1]],pretab[[k]][[2]]-minE},{k,1,dim}]]
];
ListPlot[DispTab,PlotStyle->PointSize[Medium],AxesLabel->{Momentum,Energy},AspectRatio->1]

Nump=10(*=Number of particles. The length of the ring = Nump*);
IntList=Module[{lis,inilis,dim},
lis=PadLeft[#,Nump] &/@Map[Sort,Flatten[Table[IntegerPartitions[i],{i,0,Nump*2}],1]]//Parallelize;
inilis=Table[-(Nump-1)/2+(j-1),{j,1,Nump}];
dim=Length[lis];
Table[inilis+lis[[i]],{i,1,dim}]
];
DispTab=Module[{pretab,dim,minE},
pretab=Disp[Nump,Nump,#] &/@IntList//Parallelize;
dim=Length[pretab];
minE=Min[Table[pretab[[k]][[2]],{k,1,dim}]];
Join[Table[{pretab[[k]][[1]],pretab[[k]][[2]]-minE},{k,1,dim}],Table[{-pretab[[k]][[1]],pretab[[k]][[2]]-minE},{k,1,dim}]]
];
ListPlot[DispTab,PlotStyle->PointSize[Medium],AxesLabel->{Momentum,Energy},AspectRatio->1]