%
% Script to integrate and draw and phase plots for linear autonomous 2-D systems 
% Most of the obscurity comes from the need to vectorise to make the code fast!
%
clear
clf
global a b c d
a=input(' Give a: '); % Input coefficients of linear system
b=input(' Give b: ');
c=input(' Give c: ');
d=input(' Give d: ');
tmax=input(' Give tmax: '); % Give maximum time (5 or 10 works)
hold on
t=0:0.1:tmax;
xx=-10:4:10; 
yy=xx;
[xinit,yinit]=meshgrid(xx,yy); % Set up grid of initial starting points
[m,n]=size(xinit);
mn=m*n;
yin=[reshape(xinit,mn,1);reshape(yinit,mn,1)]; % Vectorise over initial points
vecfield=rhslin(0,yin); % Get rhs directions at initial points
% Plot arrows in direction of rhs vector (gives sense of evolution)
% (See `help quiver' for usage of this MATLAB-supplied function)
quiver(yin(1:mn),yin(mn+1:2*mn),vecfield(1:mn),vecfield(mn+1:2*mn))
[tout,yout]=ode45('rhslin',t,yin); % Integrate ODEs forward in time
len=length(yin)/2;
plot(yout(:,1:len),yout(:,len+1:2*len),'r')
axis([-10 10 -10 10])
axis square
hold on
[tout,yout]=ode45('rhslin',-t,yin); % Integrate ODEs backward in time 
len=length(yin)/2;
plot(yout(:,1:len),yout(:,len+1:2*len),'r')
hold off