    N = 32;
    h = 1/N;
    
    % Option 1:
    % N-1 points per row, N-1 points per column
    % 2nd order FD [-1 4/3 -1/3] at Neumann boundary
    
    e = ones(N-1,1);
    I = eye(N-1);
    T = spdiags([e -4*e e], [-1 0 1], N-1,N-1);
    S = spdiags([e e], [-1 1], N-1,N-1);
    S(1,1) = 4/3; S(1,2) = 2/3;
    S(N-1,N-2) = 2/3; S(N-1,N-1) = 4/3;
    A = (kron(I,T) + kron(S,I))/(h^2);
    
    % Option 2:
    % N-1 points per row, N+1 points per column
    % ghost point FD [-1 0 1] at Neumann boundary
    
    f = ones(N+1,1);
    J = eye(N+1);
    T = spdiags([e -4*e e], [-1 0 1], N-1,N-1);
    S = spdiags([f f], [-1 1], N+1,N+1);
    S(1,2) = 2;
    S(N+1,N) = 2;
    B = (kron(J,T) + kron(S,I))/(h^2);
    
    
    % Full matrices
    A = full(A);
    B = full(B);
    
    % Jacobi matrices
    DA = diag(diag(A));
    GA = inv(DA)*(DA - A);
    DB = diag(diag(B));
    GB = inv(DB)*(DB - B);
   
    % Spectral radii and norms: inv(A) and inv(B)
    rho_invA = 1/min(abs(eig(A)))
    norm_invA = norm(inv(A),2)
    rho_invB = 1/min(abs(eig(B)))
    norm_invB = norm(inv(B),2)
    oneoverpisquared = 1/pi^2
    
    % Spectral radii and norms: GA and GB
    rho_GA = max(abs(eig(GA)))
    norm_GA = norm(GA,2)
    rho_GB = max(abs(eig(GB)))
    norm_GB = norm(GB,2)
    
    % Normality tests
    nonnormal_A = norm((A')*A - A*(A'),2)
    nonnormal_B = norm((B')*B - B*(B'),2)
    nonnormal_GA = norm((GA')*GA - GA*(GA'),2)
    nonnormal_GB = norm((GB')*GB - GB*(GB'),2)
    
    nonsymmetric_A = norm(A' - A,2)
    nonsymmetric_GA = norm(GA' - GA,2)

    
    
    
    
    