Sarah Adelman MATH246 Extra Credit
%Phase plane as u goes from -2 to 2
for u=-2:.1:2
A=[u+1, -u; u ,u-1];
pplane(A)
hold on
end
function pplane(A)
% pplane.m Plots the eigenvectors, eigenvalues, and representative
% phase portraits for
%
% d[x; y]/dt = A[x; y]
%
% for -1<x<1, -1<y<1. Called as
%
% pplane(A)
%
% where A is the 2x2 coefficient array. Calls lodesolver.m,
% so be sure that function is in the current path
[U,S] = eig(A);
lam = diag(S);
if lam(1) == lam(2)
% manually set u2 to be orthogonal to u1
U(:,2) = [-U(2,1); U(1,1)];
end
figure
if isreal(lam(1))
plot([-U(1,1) U(1,1)],[-U(2,1) U(2,1)],'--', ...
[-U(1,2) U(1,2)],[-U(2,2) U(2,2)],'Linewidth',2)
legend('u_1','u_2')
else
plot(0,0,'o','MarkerSize',12,'LineWidth',2)
end
xlabel('x'), ylabel('y')
grid on
trA = A(1,1)+A(2,2);
detA = A(1,1)*A(2,2) - A(1,2)*A(2,1);
title(['tr(A) = ',num2str(trA),' det(A) = ',num2str(detA)])
tmax = sqrt(abs(lam(1))^2+abs(lam(2))^2)/2;
incr = 1/max(abs(lam))/20;
t = [0:0.05:tmax];
hold on
for i = 1:40
z0 = 2*rand(2,1)-[1;1];
z = lodesolver(A,t,z0);
