using LinearAlgebra, PyPlot

A = [   0  -9   5  -4   7
       -2   2   0  -5  -9
        4   9   7   5   2
        5   0  -7   5  -4
        2   1  -3  -8   0   ]

λ = eigvals(A)

# a little function to plot the complex axes nicely
function reim_axes(ax)
    ax.axis("equal")

    # Move left y-axis and bottim x-axis to centre, passing through (0,0)
    ax.spines["left"].set_position("center")
    ax.spines["bottom"].set_position("center")

    # Eliminate upper and right axes
    ax.spines["right"].set_color("none")
    ax.spines["top"].set_color("none")

    # Show ticks in the left and lower axes only
    ax.xaxis.set_ticks_position("bottom")
    ax.yaxis.set_ticks_position("left")

    xlabel(L"\operatorname{Re} \lambda")
    ax.xaxis.set_label_coords(1.04, 0.52)

    ylabel(L"\operatorname{Im} \lambda", rotation=0)
    ax.yaxis.set_label_coords(0.5, 1.02)
end

plot(real.(λ), imag.(λ), "bo", label=L"$\lambda$ of real $A$")

reim_axes(gca())
legend()

X = eigvecs(A)
X[:,1:2] # the first two eigenvectors, corresponding to λ ≈ -5.6 ± 8.1i

B = [ -4+6im   4+3im   2+5im   6+9im   0+0im
      -3+4im   4+3im  -6-3im   5+0im   9+5im
       8+3im  -3-6im   6-3im  -9-5im   4+8im
       9+9im   9+0im   2-1im   5+2im   7+9im
      -4-9im   0+0im   6+4im  -2-6im  -4-5im ]

λ_B = eigvals(B)

plot(real.(λ_B), imag.(λ_B), "ro", label=L"$\lambda$ of complex $B$")

reim_axes(gca())
legend()

x = [1,2,3,4,5]

c = X \ x

xn = reduce(vcat, [(A^n * float(x))' for n = 1:10])

plot(xn, "o-")
ylabel(L"(A^n x)_k", size=20)
xlabel(L"matrix power $n$", size=15)

abs.(λ)

plot(xn ./ abs(λ[end]).^(1:10), "o-")
ylabel(L"(A^n x)_k / |\lambda_5|^n", size=20)
xlabel(L"matrix power $n$", size=15)

xnmore = reduce(vcat, [((A/λ[end])^n * float(x))' for n = 1:50])
plot(xnmore, "o-")
ylabel(L"(A^n x)_k / |\lambda_5|^n", size=20)
xlabel(L"matrix power $n$", size=15)