%
%  code9Fa.m -- 01 nov 02 -- djm
%
%  direction field plotting code
%

%  clear workspace & graphics ("help clear", etc)

clear;  close all

%  define matrix

disp(' ')
disp('matrix cases')
disp('  1)  decaying oscillations')
disp('  2)  growing oscillations')
disp('  3)  decaying exponentials')
disp('  4)  growing exponentials')
disp('  5)  exponential saddle')
disp('  6)  elliptical orbits')
disp('  0)  no more cases')
disp(' ')
Ca = input('case -> ');

switch Ca
    case{1}
        A = [-0.5 1.0 ; -1.0 -0.5]
    case{2}
        A = [6 -1 ; 5 4]
    case{3}
        A = [-1 1 ; -2 -4]
    case{4}
        A = [5 -1 ; 3 1]
    case{5}
        A = [0 1 ; 2 1]
    case{6}
        A = [1.5 -2.5 ; 2.5 -1.5]
    otherwise
        return
end

[ve , ei] = eig(A);
eigenvalues = diag(ei).'
first_eigenvector = ve(:,1)
second_eigenvector = ve(:,2)

%  make grid matrices

dx  = 0.2;  L = 2;
xt1 = -L-dx/2:dx:L+dx/2;
xt2 = -L-dx/2:dx:L+dx/2;

[xg1,xg2] = meshgrid(xt1,xt2);

xp1 = A(1,1)*xg1 + A(1,2)*xg2;
xp2 = A(2,1)*xg1 + A(2,2)*xg2;

r = sqrt(xp1.^2 + xp2.^2);

figure(1); clf; hold on
quiver(xt1,xt2, xp1./r, xp2./r,0.25,'k')
quiver(xt1,xt2,-xp1./r,-xp2./r,0.25,'k.')

axis equal
axis([-1 1 -1 1]*L)

%  label plot ("\bf"-prefix makes text boldface)

title(['\bf direction field'])
xlabel('\bf x_1(t)-axis')
ylabel('\bf x_2(t)-axis')