# Find parameter which solves a 9-by-9 homogeneous system of linear equations, NSolve gives some incorrect answers

I have a set of 9 linear homogeneous simultaneous equations which depend on 2 parameters, p and x. For a chosen value of p, I aim to calculate the smallest 3 values of x which satisfy the simultaneous equations. I run into trouble when plotting x as a function of p, which I think is down me using NSolve incorrectly.

## Background functions for the problem

The elements which form the coefficient matrix for this problem are in turn calculated from the eigenfunctions/values from another square matrix, given below:

ClearAll[eigs, value1, func1, value2, func2]
g = 1.006;
ham[p_] := SparseArray[{Band[{1, 1}] -> Table[i^2, {i, -9, 9}], Band[{2, 1}] -> p*100*(1 - g), Band[{1, 2}] -> p*100*(1 - g), Band[{3, 1}] -> p*(g/2), Band[{1, 3}] -> p*(g/2)}, {19, 19}]
(* 19x19 Matrix from which I find eigenfunctions and eigenvalues which will later contribute to my system of simultaneous equations. I will need eigenfunctions and values for both ham[p] and ham[-p]. *)

eigs[p_] := eigs[p] = SortBy[Transpose[Eigensystem[ham[p]]], First] (* Eigenvalues of the above matrix, ordered with smallest eigenvalue listed first. *)
value1[p_, n_] := value1[p, n] = eigs[p][[n]][[1]] (* Eigenvalues *)
func1[p_, n_] := func1[p, n] = Chop[(1/Sqrt[2*Pi])*eigs[p][[n]][[2]]] . Table[Exp[I*j*\[CurlyPhi]], {j, -9, 9}] (* Function constructed as an exponential series with the eigenvector elements as coefficients. *)
value2[p_, n_] := value2[p, n] = value1[-p, n] (* Eigenvalues for matrix ham[-p] *)
func2[p_, n_] := func2[p, n] = func1[-p, n] (* Similar to func1 but for eigenvectors from ham[-p] *)


## Coefficient matrix for system of linear equations

I define here the coefficient matrix for the homogeneous set of linear equations I wish to solve. Its a 9-by-9 matrix, which I input manually. Each element is formed by a numerical integral involving funk1 & func2, multiplied by a prefactor involving value1 and value2. I apologise in advance if this is an eyesore:

Clear[coeff]

me[(p_)?NumericQ, x_, (n0_)?NumericQ, (M0_)?NumericQ, (n_)?NumericQ,(M_)?NumericQ] := Chop[
((0.4*Pi)/(value1[p, n] + value2[p, M] - x))*
NIntegrate[func1[p, n]*func2[p, M]*Conjugate[func1[p, n0]]*Conjugate[func2[p, M0]], {\[CurlyPhi], -Pi, Pi},
AccuracyGoal -> 6, Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule", "SymbolicProcessing" -> 0}]];

coeff[p_, x_] := coeff[p, x] = {
{me[p, x, 1, 1, 1, 1] - 1, me[p, x, 1, 1, 1, 2], me[p, x, 1, 1, 1, 3], me[p, x, 1, 1, 2, 1], me[p, x, 1, 1, 2, 2], me[p, x, 1, 1, 2, 3], me[p, x, 1, 1, 3, 1], me[p, x, 1, 1, 3, 2], me[p, x, 1, 1, 3, 3]},
{me[p, x, 1, 2, 1, 1], me[p, x, 1, 2, 1, 2] - 1, me[p, x, 1, 2, 1, 3], me[p, x, 1, 2, 2, 1], me[p, x, 1, 2, 2, 2], me[p, x, 1, 2, 2, 3], me[p, x, 1, 2, 3, 1], me[p, x, 1, 2, 3, 2], me[p, x, 1, 2, 3, 3]},
{me[p, x, 1, 3, 1, 1], me[p, x, 1, 3, 1, 2], me[p, x, 1, 3, 1, 3] - 1, me[p, x, 1, 3, 2, 1], me[p, x, 1, 3, 2, 2], me[p, x, 1, 3, 2, 3], me[p, x, 1, 3, 3, 1], me[p, x, 1, 3, 3, 2], me[p, x, 1, 3, 3, 3]},
{me[p, x, 2, 1, 1, 1], me[p, x, 2, 1, 1, 2], me[p, x, 2, 1, 1, 3], me[p, x, 2, 1, 2, 1] - 1, me[p, x, 2, 1, 2, 2], me[p, x, 2, 1, 2, 3], me[p, x, 2, 1, 3, 1], me[p, x, 2, 1, 3, 2], me[p, x, 2, 1, 3, 3]},
{me[p, x, 2, 2, 1, 1], me[p, x, 2, 2, 1, 2], me[p, x, 2, 2, 1, 3], me[p, x, 2, 2, 2, 1], me[p, x, 2, 2, 2, 2] - 1, me[p, x, 2, 2, 2, 3], me[p, x, 2, 2, 3, 1], me[p, x, 2, 2, 3, 2], me[p, x, 2, 2, 3, 3]},
{me[p, x, 2, 3, 1, 1], me[p, x, 2, 3, 1, 2], me[p, x, 2, 3, 1, 3], me[p, x, 2, 3, 2, 1], me[p, x, 2, 3, 2, 2], me[p, x, 2, 3, 2, 3] - 1, me[p, x, 2, 3, 3, 1], me[p, x, 2, 3, 3, 2], me[p, x, 2, 3, 3, 3]},
{me[p, x, 3, 1, 1, 1], me[p, x, 3, 1, 1, 2], me[p, x, 3, 1, 1, 3], me[p, x, 3, 1, 2, 1], me[p, x, 3, 1, 2, 2], me[p, x, 3, 1, 2, 3], me[p, x, 3, 1, 3, 1] - 1, me[p, x, 3, 1, 3, 2], me[p, x, 3, 1, 3, 3]},
{me[p, x, 3, 2, 1, 1], me[p, x, 3, 2, 1, 2], me[p, x, 3, 2, 1, 3], me[p, x, 3, 2, 2, 1], me[p, x, 3, 2, 2, 2], me[p, x, 3, 2, 2, 3], me[p, x, 3, 2, 3, 1], me[p, x, 3, 2, 3, 2] - 1, me[p, x, 3, 2, 3, 3]},
{me[p, x, 3, 3, 1, 1], me[p, x, 3, 3, 1, 2], me[p, x, 3, 3, 1, 3], me[p, x, 3, 3, 2, 1], me[p, x, 3, 3, 2, 2], me[p, x, 3, 3, 2, 3], me[p, x, 3, 3, 3, 1], me[p, x, 3, 3, 3, 2], me[p, x, 3, 3, 3, 3] - 1}
};


## Using NSolve to find x values for a given value of p

I finally get to the crux of my problem. I numerically solve by setting the determinant of the coefficient matrix to 0, and restrict x to between -10 and 10 (I have reasonable certainty that my solutions will not be outside this range for the range of p's I intend to work with).

xsolution[p_] := Sort[ NSolve[Det[coeff[p, x]] == 0. && (-10 <= x <= 10), x] ] (* Numerically solve which values of x yields Det[coeff]=0, and sorts the values. *)


Plotting what I believe are the 3 lowest values of x that I've found (as a function of p), I find that I have some odd-looking discontinuities, especially in the upper curve.

Plot[{
xsolution[p][[1]][[1]][[2]],
xsolution[p][[2]][[1]][[2]],
xsolution[p][[3]][[1]][[2]]
}, {p, -2, 2}, MaxRecursion -> 4, PlotPoints -> 5]


What's more, I do not trust my results from NSolve! For example, the 2nd lowest value of x below (for p=1) does not yeild Det[coeff]=0 as expected. What am I doing wrong?!

xsolution[1]


yields output:

{{x -> -1.321793518956363}, {x -> -1.2613922782902645}, {x -> -0.5449889067896768}, {x -> -0.4556837316778105}, {x -> -0.455683729053839}, {x -> 0.34035978790527527}, {x -> 0.48514450207601595}, {x -> 0.682370748049176}, {x -> 0.6824136788022568}, {x -> 0.7017597948788745}, {x -> 0.7564607365305215}, {x -> 0.7565115382635393}, {x -> 1.3495866699741845}, {x -> 1.4698703476512056}, {x -> 1.4881232169456888}, {x -> 1.4881232372192312}, {x -> 1.690848785958313}, {x -> 1.69087512852734}, {x -> 2.5257845208148724}, {x -> 2.700249525479737}, {x -> 2.7003184957110764}, {x -> 2.7062611577155558}, {x -> 2.9030703703488534}, {x -> 2.903070397152075}, {x -> 2.9030704239717173}}


but

Det[coeff[1, -1.2613922782902645]]


gives

-28019.791238792383


Any advice would be greatly appreciated! I'm relatively inexperienced in Mathematica (as probably abundantly clear in my code above), so any help with my problems - or advice on speeding up these clunky calculations - would be very much appreciated!! Many thanks.

• Welcome to Mma.SE. Thanks for taking the tour. Be sure you have learned about asking and what's on-topic. Always edit if improvable. Showing due diligence, giving brief context, including minimal working examples of code and data in formatted form and so on help us to help you and in most cases inspires great answers. The site depends on participation, as you receive give back: vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. – rhermans Aug 21 '18 at 16:43