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

%Part 1 - Direction field


P =[1 -1;1 3];
[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,1,1000);
y1=[exp(2*t);-exp(2*t)];
y2=[(t+1).*exp(2*t);-(t+2).*exp(2*t)];
c1=3; c2=-2;
y=c1*y1+c2*y2;
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(-3,3,1000);
y1=[exp(2*t);-exp(2*t)];
y2=[(t+1).*exp(2*t);-(t+2).*exp(2*t)];

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

C1=(-2:.5:2); C2=(-2:.5:2);
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([-3 3 -3 3])
box


