%Math 316 - Spring 2017
%MATLAB code for System2 of Chapter 7 
%Written by M. Motamed @ UNM 
%===========================================================
%Direction field and phase portrait 
%for the ODE system: y'=Py with P =[-0.5 1;-1 -0.5]
%IC: y0=[1;1]
%===========================================================

%Part 1 - Direction field


P=[-.5 1;-1 -.5];
[Y1,Y2] = meshgrid(-3:.25:3,-3:.25:3);
Y1=Y1(:); Y2=Y2(:);

N=length(Y1);

DY1=[]; DY2=[];
for i=1:N 
    y=[Y1(i);Y2(i)];
    dy=P*y;
    DY1(i,1)=dy(1); DY2(i,1)=dy(2);
    
end

L=sqrt(DY1.^2 + DY2.^2);

figure
set(gca,'fontsize', 20); 
quiver(Y1,Y2,DY1./L,DY2./L)

xlabel('y_1')
ylabel('y_2')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 2 - Plotting the solution
t=linspace(0,10,1000);
u=[cos(t).*exp(-0.5*t);-sin(t).*exp(-0.5*t)];
v=[sin(t).*exp(-0.5*t);cos(t).*exp(-0.5*t)];
c1=1; c2=1;
y=c1*u+c2*v;
figure
set(gca,'fontsize', 20); 
plot(t,y(1,:),'b-',t,y(2,:),'r-')
legend('y1(t)','y2(t)')
xlabel('t')
ylabel('y(t)')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 3 - Phase portrait
t=linspace(-10,10,1000);
u=[cos(t).*exp(-0.5*t);-sin(t).*exp(-0.5*t)];
v=[sin(t).*exp(-0.5*t);cos(t).*exp(-0.5*t)];

figure
set(gca,'fontsize', 20);
hold on

C1=(-2:0.5:2); C2=(-2:.5:2);
for i=1:length(C1)
    for j=1:length(C2)
        c1=C1(i); c2=C2(j);
        y=c1*u+c2*v;
        plot(y(1,:),y(2,:),'k-')
    end
end
xlabel('y_1')
ylabel('y_2')
axis equal
axis([-3 3 -3 3])




