The spectrum of tridiagonal Toeplitz matrices is known analytically (Tridiagonal-Toeplitz-Matrix – Wikipedia) and is (for non-zero off-diagonals) simple and real. Julia’s eigen function returns for small non-symmetric Toeplitz matrices the correct real spectrum, but for larger sizes fails to do so. E.g.

n = 200
p = 0.7
q = 0.4
A = Tridiagonal([q for _ in 1:(n-1)], zeros(n), [p for _ in 1:(n-1)])
A = Array(A)
maximum(imag.(eigvals(A)))

Shows non-zero imaginary parts significantly larger than machine precision. The same happens in numpy. Changing to BigFloat helped for sizes of n around 200. I did not test for larger n, but would guess that at sufficiently large n it will fail as well.

The eigenvalues of these matrices are super ill conditioned. I spent some time on these matrices with @sverek a while ago. Sven-Erik knows a lot more about these than I do.

The condition number of the eigenvalue is roughly reciprocal of the inner product between the right and the left (normalized) eigenvectors for that value. (I think I have the condition number from Nick HIgham’s book). So you have

julia> n = 20;
julia> A = Tridiagonal([q for _ in 1:(n-1)], zeros(n), [p for _ in 1:(n-1)]);
julia> F = eigen(Array(A));
julia> vr = (A - F.values[end]*I)\randn(n) |> t -> t/norm(t);
julia> vl = (A' - F.values[end]*I)\randn(n) |> t -> t/norm(t);
julia> abs(inv(vl'vr))
6.737730466564851
julia> n = 200;
julia> A = Tridiagonal([q for _ in 1:(n-1)], zeros(n), [p for _ in 1:(n-1)]);
julia> F = eigen(Array(A));
julia> vr = (A - F.values[end]*I)\randn(n) |> t -> t/norm(t);
julia> vl = (A' - F.values[end]*I)\randn(n) |> t -> t/norm(t);
julia> abs(inv(vl'vr))
7.360334831197046e19

so n doen’t have to get large before small perturbations will be heavily amplified in the eigenvalues.

I can send you some current references regarding this (and some ways how to deal with it). By the way, the spectrum does not have to be real (e.g., just set q=-0.4).

using LinearAlgebra
using Plots
using Pseudospectra
n = 1000
p = 0.7
q = 0.4
A = Tridiagonal([q for _ in 1:(n-1)], zeros(n), [p for _ in 1:(n-1)])
A = Array(A)
λ = eigvals(A)
spectralportrait(A)

Thank you for all your great answers. I was not aware that finding the eigenvalues of tridiagonal Toeplitz matrices is such an ill-conditioned problem. I will mark @andreasnoack post as answer, but I am very grateful for all your comments, especially @zdenek_hurak pointing to the pseudo spectrum.

@zdenek_hurak@rafael.guerra I did my tests in Windows 10. On WSL 2, Ubuntu 20.04 I get non-zero imaginary parts for n=200 as well. My specs are for
Windows

Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
JULIA_NUM_THREADS = 4

and WSL 2 Ubuntu:

julia> versioninfo()
Julia Version 1.6.3
Commit ae8452a9e0 (2021-09-23 17:34 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, skylake)