t=linspace(0,8,200);
A0=[[-1 0]; [0 -1]]
A1=[[-1 1]; [0 -1]]
A2=[[-1 5] ; [0 -2]]
for j=1:length(t)
    f0(j)=norm(expm(t(j)*A0));
    f1(j)=norm(expm(t(j)*A1));
   f2(j)=norm(expm(t(j)*A2));
end
figure, plot(t,f0,'r',t,f1,'g',t,f2,'b')
title('Matrix exponentials, e^{tA}')
pause
legend('A0', 'A1','A2')
pause

dx=.000001
figure
t=linspace(0,2*pi,400);
%Perturbation of the n x n Jordan block by dx in the (n,1) position
for j=1:20
    n=2*j;
    A=gallery('forsythe',n,dx,0);
    da=eig(A);
    scatter(real(da),imag(da))
    hold on
    r=dx^(1/n);
    plot(r*cos(t),r*sin(t))
    axis equal
    hold on
end
title(['Perturbations in the (n,1) spot of the n x n Jordan Block, n = 2 to ', num2str(n)])
hold off
pause
m=50
figure
%Random perturbation of the m x m Jordan block
for l=1:50
    B=dx*rand(m)+gallery('forsythe',m,0,0);
    db=eig(B);
    scatter(real(db),imag(db))
    axis equal
    title(['Random Perturbations of the ',num2str(m),' x ',num2str(m),' Jordan Block'])
    hold on
end
hold off
pause
p=40;
dy=.0005
%Perturbation of a tridiagonal matrix with 1 below the diagonal and -od
%above. We can see the transition from non-normal toward normal behavior.
%Note that A=gallery('tridiag',p,1,0,-1) is a skew symmetric matrix.
%Hence iA is self adjoint and therefore A is normal.
eps=.025;
figure
for l=1:41
    od=eps*(l-1);
    B=dy*rand(p)+gallery('tridiag',p,1,0,-od);
    db=eig(B);
    scatter(real(db),imag(db))
    dd=eig(zeros(p)+gallery('tridiag',p,1,0,-od));
    hold on
    scatter(real(dd),imag(dd),'x','r')
    xlim([-2.5,2.5]) 
    ylim([-2.5,2.5])
    title(['Spectrum of perturbations of (1,0,',num2str(od),') tridiagonal matrix, dy = ',num2str(dy)])
    hold off
    legend('Perturbed spectrum','Unperturbed spectrum')
    pause
end
%In the previous example there seemed to be something fishy about the
%"exact" eigenvalues when the upper diagonal was very close zero. We examine that
%here.
pause
eps = .025
figure
for l=1:40
    od=eps;
    dd=eig(zeros(2*l)+gallery('tridiag',2*l ,1,0,-od));
    scatter(real(dd),imag(dd),'x','r')
    hold on
    scatter(zeros(1,2*l),2*sqrt(eps)*cos((1:(2*l))*pi/(2*l+1)),'g');
    xlim([-0.5,0.5]) 
    ylim([-0.5,0.5])
    title([num2str(2*l),' x ',num2str(2*l),'(1,0,',num2str(-eps),') tri-diag matrix'])
    hold off
    pause
end
