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

%Part 1 - Direction field

P=[1 1;4 1];
[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
hold on
set(gca,'fontsize', 20); 
quiver(Y1,Y2,DY1./L,DY2./L)

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

plot([-4 4],[0 0],'r-')
plot([0 0],[-4 4],'r-')


%if 1==2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Part 2 - Phase portrait

t=linspace(-5,5,10000);
y1=[exp(t);-2*exp(t)];
y2=[exp(3*t);2*exp(3*t)];

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

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

%end