function P1()

clc
n1 = floor(30 + 10*rand());
n2 = floor(30 + 10*rand());
a = -(5*(1+rand()));
b = -(5*(1+rand()));
c = -4*(a+b);
A = MakeA(n1,n2,a,b,c);
eigsA = eig(A);
TrueMin = min(eigsA);
TrueMax = max(eigsA);
[LambdaMin,LambdaMax] = EigsA(n1,n2,a,b,c);
fprintf('n1 = %2d n2 = %2d  a = %5.3f  b = %5.3f  c = %5.3f\n\n',n1,n2,a,b,c)
fprintf('TrueMin = %15.6f\n',TrueMin)
fprintf('MyMin   = %15.6f\n',LambdaMin)
fprintf('TrueMax = %15.6f\n',TrueMax)
fprintf('MyMax   = %15.6f\n',LambdaMax)
MyRho = SpectralRad(n1,n2,a,b,c);
TrueRho = max(abs(eig(A-c*eye(n1*n2,n1*n2))))/c;
fprintf('TrueRho = %15.6f\n',TrueRho)
fprintf('MyRho   = %15.6f\n',MyRho)

function A = MakeA(n1,n2,a,b,c)
% Explicit formation of the A-Matrix
In1 = eye(n1,n1);
In2 = eye(n2,n2);
En1 = diag(ones(n1-1,1),-1) + diag(ones(n1-1,1),1);
En2 = diag(ones(n2-1,1),-1) + diag(ones(n2-1,1),1);
A =  kron(In1,c*In2+b*En2) + kron(En1,b*In2+a*En2);

function [LambdaMin,LambdaMax] = EigsA(n1,n2,a,b,c)
% The smallest and largest eigenvalues of the A-matrix
lambda = zeros(n1,n2);
for i=1:n1
    for j=1:n2
        alfa = cos(j*pi/(n2+1));
        beta = cos(i*pi/(n1+1));
        lambda(i,j) = c + 2*b*(alfa+beta) + 4*a*alfa*beta;
        lambda(i,j) = -2*b*((1-alfa) + (1-beta)) - 4*a*(1-alfa*beta);
    end
end
% Assume that b and c are negative.
% Smallest eigenvalue is lambda(1,1), for then alfa and beta are maximized.
% By making alfa and beta as negative as possible, we get the largest
% eigenvalue, i.e., lambdaIn1,n2).
LambdaMin = lambda(1,1);
LambdaMax = max(lambda(n1,1),lambda(1,n2));

function rho = SpectralRad(n1,n2,a,b,c)
% The A-matrix has diagonal equal to cI, so the Jacobi iteration matrix is
%      -(kron(In1,b*En2) + kron(En1,b*In2+a*En2))/c
% Lets compute its eigs
lambda = zeros(n1,n2);
for i=1:n1
    for j=1:n2
        alfa = cos(j*pi/(n2+1));
        beta = cos(i*pi/(n1+1));
        lambda(i,j) = -(2*b*(alfa+beta) + 4*a*alfa*beta)/c;
    end
end
rho = max(abs(lambda(1,1)));
 
    