% Uses quiver command to plot directions as field of arrows, for a
% two-dimensional state space with coordinates x1,x2

% !!! TYPE help quiver to get more information on how this command is used

% Uses my_odefun.m 
% This must be a separate function (stored in same directory) with the
% standard ODE syntax:

% function F=my_odefun(t,x).  
    %x is the state vector
    %F is the velocity vector
    %t is the time.  Note, you have to use this t argument, even if it is
    %not used in defining the velocity!

    
t=0;  %MUST define a value of time, even if it is note used in defining the velocity    
    
%define range of x1 and x2 values over which make the plot
x1min=-1;
x1max=1;
x1spacing=0.2;
x1list=x1min:x1spacing:x1max ;

x2min=-1;
x2max=1;
x2spacing=0.2;
x2list=x1min:x1spacing:x1max ;


%OK, we are going to define a grid of x1 and x2 values, based on the lists
%above.  The (i,j) point on that grid will correspond to elements (i,j) in
%four matrices that we are about to make.

%x1 values increase across ROWS (horizontal direction, index j)
%x2 values increase down COLUMNS (vertical direction, index i)

[x1matrix,x2matrix]=meshgrid(x1list,x2list);

%That makes two matrices.  Element (i,j) of x1matrix is the x1 coordinate
%for that (i,j). Element (i,j) of x2matrix is the x2 coordinate
%for that (i,j).  


%make two more "dummy" matrices
dx1dt_matrix=zeros(size(x1matrix));
dx2dt_matrix=zeros(size(x2matrix));


for i=1:length(x2list)
    for j=1:length(x1list)
        
    %here, my_odefun will return a vector F of velocities in x1 and x2 directions    
    F=my_odefun(0,[x1matrix(i,j);x2matrix(i,j)]);

    %I fill in these two matrices with the corresponding velocities
    dx1dt_matrix(i,j)=F(1);
    dx2dt_matrix(i,j)=F(2);
    end
end


%finally, the quiver command makes a plot of arrows according to these
%velocities.

figure
set(gca,'FontSize',16)
quiver(x1matrix,x2matrix,dx1dt_matrix,dx2dt_matrix,.5,'LineWidth',1)
xlabel('x1','FontSize',20);ylabel('x2','FontSize',20)

pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Next, in rescaled variables

vlist=0:.1:1;
slist=0:.1:1;

[vmatrix,smatrix]=meshgrid(vlist,slist);

dvdt_matrix=zeros(size(vmatrix));
dsdt_matrix=zeros(size(vmatrix));

for i=1:length(slist)
    for j=1:length(vlist)
    dxdt=michaelis_mentin_rescaled_odefun(0,[vmatrix(i,j);smatrix(i,j)]);
    dvdt_matrix(i,j)=dxdt(1);
    dsdt_matrix(i,j)=dxdt(2);
    end
end


figure
set(gca,'FontSize',16)
quiver(vmatrix,smatrix,dvdt_matrix,dsdt_matrix,'LineWidth',1)
xlabel('v','FontSize',20);ylabel('s','FontSize',20)