This is the matrix approximating the Laplacian plus the boundary conditions u(0) = u(π) = 0, with n = 5 subdivisions of the interval [ 0, π].  We know that the true eigenvalues are 
(Subscript[λ, k]) = -(k^2) and the eigenfunctions are (Subscript[u, k])= sin(k π x/π) = sin(k x), for k = 1, 2, 3, ..., n

Lunscaled =  {{-3,1,0,0, 0,0,0,0},{1,-2,1,0, 0,0,0,0},{0,1,-2,1, 0,0,0,0}, {0, 0,1,-2,1,0,0,0}, {0,0, 0,1,-2,1,0,0},{0,0,0, 0,1,-2,1,0},{0, 0,0,0,0,1,-2,1},{0,0,0,0, 0,0,1,-3}};
Print["L = (64/(π^2)) ", MatrixForm[Lunscaled]]
L= 64./Pi^2 Lunscaled;  (* I put the decimal point in 64. to make it a floating point matrix, not a symbolic matrix *)

Eigenvalues[L]  (* This lists them from largest in absolute value to smallest *)
Reverse[Eigenvalues[L]]  (* List them the reverse order *)

Now, consider the eigenvalues and eigenvectors of the n x n matrix

n =5; (* make this at least 5.  For n < 21 the L matrix is displayed *)
L = Table[0, {n}, {n}];
L[[1,1]] = -3; (* -2 + u0/u1 in general *)
L[[1, 2]] = 1;
L[[n, n-1]] = 1;
L[[n,n]] = -3; (* -2 + un+1/un *)
Do[L[[i, i]] = -2;
L[[i, i-1]] = 1;
L[[i, i+1]] = 1,
{i, 2, n-1}];
If[n < 21, Print["Matrix approximating the Laplacian with BCs u(0) = u(π) = 0:\nL = ", ((n^2)/(π^2)),MatrixForm[L] ]];
L = n^2./Pi^2 L;  (* the 2. makes it a floating point matrix. Try taking it out when n is small *)
evals = Reverse[Eigenvalues[L]];
Print["The first 5 eigenvalues of L with n = ", n, " are\n", Take[evals,5] ];
evecs = Reverse[Eigenvectors[L]];
pd = ListPlot[evecs[[2]], Joined -> True];
a = evecs[[2,1]]/Sin[.5* 2 Pi/n];
ps = Plot[a Sin[(i-.5) 2 Pi/n], {i, .5, n + .5}, PlotStyle -> Red, Ticks ->{{{(n+1)/2,π/2}, {n+1/2,π}}, Automatic}] ;
Show[ps, pd, PlotLabel -> "Second Eigenfunction"]

Now, find the matrix for second derivative with Boundary conditions u’(0) = u’(π) = 0.

n =5; (* make this at least 5.  It takes about a second for n = 500 *)
L = Table[0, {n}, {n}];
L[[1,1]] = -1; (* -2 + u0/u1 in general *)
L[[1, 2]] = 1;
L[[n, n-1]] = 1;
L[[n,n]] = -1; (* -2 + un+1/un *)
Do[L[[i, i]] = -2;
L[[i, i-1]] = 1;
L[[i, i+1]] = 1,
{i, 2, n-1}];
If[n < 21, Print["Matrix approximating the Laplacian with Neumann Boundary conditions:\nL = ", ((n^2)/(π^2)),MatrixForm[L] ]];
L = n^2/Pi^2 L;
evals = Reverse[Eigenvalues[L//N]];
Print["The first 5 eigenvalues of L with n = ", n, " are\n", Take[evals,5] ];
evecs = Reverse[Eigenvectors[L//N]];
pd = ListPlot[evecs[[2]], Joined -> True];
a = evecs[[2,1]]/Cos[.5* Pi/n];
ps = Plot[a Cos[(i-.5)  Pi/n], {i, .5, n + .5}, PlotStyle -> Red, Ticks ->{{{(n+1)/2,π/2}, {n+1/2,π}}, Automatic}] ;
Show[ps, pd, PlotLabel -> "Second Eigenfunction"]

Now, do the matrix for the BC u(0) = 0, u(π) + u’(π) = 0.
First, compute the second eigenvalue numerically.  Setting eval = k^2,  we get u(x) = sin(k x) and we we need to solve a “transcendental” equation for k to satisfy the BC at x = π.

Plot[Sin[k π] + k Cos[k π], {k, 0, 3}]

FindRoot[Sin[k π] + k Cos[k π] == 0, {k, 0.8}]  (* This returns an "assignment rule", not the actual value of the root *)
k1 = k/.FindRoot[Sin[k π] + k Cos[k π] == 0, {k, 0.8}]
k2 = k/.FindRoot[Sin[k π] + k Cos[k π] == 0, {k, 1.7}]
k3 = k/.FindRoot[Sin[k π] + k Cos[k π] == 0, {k, 2.6}]

To get the bottom right element of the L matrix, we need to solve for (Subscript[u, n+1])/(Subscript[u, n]) with the assumption that u + u’ = 0 at the boundary.  I’ll do this at the board and get ((Subscript[u, n+1])/(Subscript[u, n]))= ((2n-π)/(2n+π)).

n =5; (* make this at least 5.  It takes about a second for n = 500 *)
L = Table[0, {n}, {n}];
L[[1,1]] = -3 ;(* -2 + u0/u1 in general *)
L[[1, 2]] = 1;
L[[n, n-1]] = 1;
L[[n,n]] = -2 + (2n-π)/(2n+π); (* -2 + un+1/un *)
Do[L[[i, i]] = -2;
L[[i, i-1]] = 1;
L[[i, i+1]] = 1,
{i, 2, n-1}];
If[n < 21, Print["Matrix approximating the Laplacian with Neumann Boundary conditions:\nL = ", ((n^2)/(π^2)),MatrixForm[L] ]];
L = n^2/Pi^2 L;
evals = Reverse[Eigenvalues[L//N]];
Print["The first 5 eigenvalues of L with n = ", n, " are\n", Take[evals,5] ];
evecs = Reverse[Eigenvectors[L//N]];
pd = ListPlot[evecs[[2]], Joined -> True];
a = evecs[[2,1]]/Sin[k2*.5* Pi/n];
ps = Plot[a Sin[k2(i-.5)  Pi/n], {i, .5, n + .5}, PlotStyle -> Red, Ticks ->{{{(n+1)/2,π/2}, {n+1/2,π}}, Automatic}] ;
Print["Second eigenvalue is -k2^2 = -", k2^2];
Show[ps, pd, PlotLabel -> "Second Eigenfunction"]

Now do the  matrix with periodic boundary conditions.

n =9; (* make this at least 9.  It takes about a second for n = 500 *)
L = Table[0, {n}, {n}];
L[[1,1]] = -2;
L[[1, 2]] = 1;
L[[1, n]] = 1;
L[[n, n-1]] = 1;
L[[n,n]] = -2;
L[[n, 1]] = 1;
Do[L[[i, i]] = -2;
L[[i, i-1]] = 1;
L[[i, i+1]] = 1,
{i, 2, n-1}];
If[n < 21, Print["Matrix approximating the Laplacian with Periodic Boundary conditions:\nL = ", ((n^2)/(((2π))^2)),MatrixForm[L] ]];
L = n^2/(2Pi)^2 L;
evals = Reverse[Eigenvalues[L//N]];
Print["The first 5 eigenvalues of L with n = ", n, " are\n", Take[evals,9] ];
evecs = Reverse[Eigenvectors[L//N]];
Print["Plot of the first 3 eigenfunctions"];
p1 = ListPlot[{evecs[[1]], evecs[[2]], evecs[[3]]}, Joined -> True]