%
% 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
global xmin xmax ymin ymax
global fswitch arcswitch
global epsilon
epsilon=input(' enter epsilon: ');
arcswitch=input(' independent variable, enter 1 for arc-length, 0 for time: ');

if(arcswitch==1)
  smax=input(' Give smax (e.g. 5): '); % Give maximum arc-length
else
  smax=input(' Give tmax (e.g. 5): ');
end

xmax=input(' Give xmax (e.g. 5): ');
fswitch=input(' enter ref. no. for ODE (tutorial 9, first three, 1,2 or 3): ');
hold on
nsq=input(' input approx. no. of trajectories (e.g. 50): ');
n=floor(sqrt(nsq));
s=0:smax/1000:smax;
xmin=-xmax;
ymin=xmin;
ymax=xmax;
delx=(xmax-xmin)/(n-1);
dely=(ymax-ymin)/(n-1);
[xinit,yinit]=meshgrid(xmin:delx:xmax,ymin:dely:ymax);
% add randomness for better coverage of the plane
xinit=xinit+0.3*(rand(n,n)-1)*(xmax-xmin); 
yinit=yinit+0.3*(rand(n,n)-1)*(ymax-ymin);

[m,n]=size(xinit);
mn=m*n;
xv=reshape(xinit,mn,1);
yv=reshape(yinit,mn,1);
yin=[xv;yv];                            % Vectorise over initial points
vecfield=rhsarc(0,yin);                 % Get rhs directions at initial points
% Plot arrows in direction of rhs vector (gives sense of evolution)
ang=atan2(vecfield(mn+1:2*mn),vecfield(1:mn));
arrowhead(xv,yv,ang,xmax/20,'b');
hold on

[sout,yout]=ode45('rhsarc',s,yin);     % Integrate ODEs forward in arc-length
len=length(yin)/2;
plot(yout(:,1:len),yout(:,len+1:2*len),'b')
set(gca,'fontsize',14)
xlabel('x')
ylabel('y')
switch fswitch
  case 1
    title('dxdt=x-y; dydt=x+y-2xy')
  case 2 
    title('dxdt=y e^y; dydt=1-x^2')
  case 3
    title('dxdt=1-xy; dydt=(x-1)y')
  otherwise
end

box
axis([xmin xmax ymin ymax])
axis square 
hold on
[sout,yout]=ode45('rhsarc',-s,yin);    % Integrate ODEs backward in arc-length 
len=length(yin)/2;
plot(yout(:,1:len),yout(:,len+1:2*len),'b')

% add additional trajectories
fprintf('add trajectory? (click outside graph box  to finish)\n');
[xtraj,ytraj]=ginput(1);
while xtraj > xmin & xtraj < xmax & ytraj > ymin & ytraj < ymax
  yin=[xtraj;ytraj];
  vecfield=rhsarc(0,yin);
  ang=atan2(vecfield(2),vecfield(1));
  [sout,youta]=ode45('rhsarc', s,yin);  % Integrate ODEs forward in arc-length
  axis([xmin xmax ymin ymax])
  hold on
  [sout,youtb]=ode45('rhsarc',-s,yin);  % Integrate ODEs backward in arc-length

  plot(youta(:,1),youta(:,2),'r')        % plot most recent trajectory in red
  plot(youtb(:,1),youtb(:,2),'r')
  arrowhead(yin(1),yin(2),ang,xmax/20,'r')  
  [xtraj,ytraj]=ginput(1);               % starting point for next trajectory
  plot(youta(:,1),youta(:,2),'b')        % plot again in blue when moving to next trajectory 
  plot(youtb(:,1),youtb(:,2),'b')
  arrowhead(yin(1),yin(2),ang,xmax/20,'b')  
end
hold off