Numerical value of Determinant far from what it is supposed to be

As you suspected when you mentioned SetPrecision, you are encountering numerical errors, probably catastrophic loss of precision when calculating the determinant; your calculations do in fact need to be carried out at higher precision.

If possible, you would want to use exact numbers in your matrix, or take advantage of the arbitrary-precision capabilities of Mathematica. For instance, we can convert all machine-precision numbers to arbitrary-precision ones with a number of digits of precision equal to that of common machine-precision numbers on your machine using SetPrecision (see also $MachinePrecision in the documentation):

det = Det[SetPrecision[mat, $MachinePrecision]];
sol = NSolve[det == 0, h];
det /. sol // PossibleZeroQ

(* Out: 
 {True, True, True, True, True, True, True, True, True, True, True,
  True, True, True, True, True, True, True, True, True, True, True,
  True, True} 
*)

As you can see, all those values of $h$ do bring your determinant reasonably close to zero, within machine-precision approximations.

As already noted, this is a polynomial eigenproblem. First, let's use SeriesCoefficient[] to extract the coefficient matrices of your $h$-matrix:

coeffs = Table[SeriesCoefficient[mat, {h, 0, k}], {k, Max[Exponent[mat, h]], 0, -1}];

You can then use PolynomialEigenvalues[] to find the eigenvalues of your $h$-matrix:

eigs = PolynomialEigenvalues[coeffs];

To check the eigenvalues returned:

ListLinePlot[Log10[Abs[Det[mat /. h -> #] & /@ eigs]], 
             Axes -> None, Frame -> True, PlotRange -> All]

and we see that the small eigenvalues are computed accurately, but the large eigenvalues are not. (Why this is so, I still have to do some research on.)

So, as a sanity check, let us compute the eigenvalues a little differently:

eigs2 = Reverse[1/PolynomialEigenvalues[Reverse[coeffs]]];
ListLinePlot[Log10[Abs[Det[mat /. h -> #] & /@ eigs2]],
             Axes -> None, Frame -> True, PlotRange -> All]