%% Perturb the diagonal element of a Hermitian matrix

M = [1 2 3 4;
    0 5 6 5;
    0 0 4 3;
    0 0 0 2];
M = (M + M')/2;

dP = diag([0 1 0 0]);

P = M+10*dP;

disp(eig(M)');
disp(eig(P)');

%% Plot the eigenvalue as functions of perturbation strength

ee = 0.001;

arange = linspace(0,pi/2-ee);
% arange = linspace(-pi/2+ee,pi/2-ee, 201);
% arange = [arange, linspace(pi/2+ee,3*pi/2-ee, 201)];

eigs = zeros(size(M,1), length(arange));

for j = 1:length(arange),
    P = M + tan(arange(j))*dP;
    eigs(:,j) = eig(P);
end;

plot(arange, eigs, '.');

axis([min(arange), max(arange), -15, 15]);
grid on;

%% Pert strength -> infinity

P = M + 100*dP;
sp = eig(P)';

disp(sp(1:3));

Pinf = M([1 3 4], [1 3 4]);
disp(eig(Pinf)');

% or define Pinf more generally
R = null( dP );
Pinf = R' * M * R;
disp(eig(Pinf)');
