% Tht1 - study the system a little bit
% plot the direction - field 
% Etc...
% 
% 


function tht2a(u)



x = u(1:2);

if( length(u) ==3)
  a(1) = u(3);
  a(2) = 14;
  a(3) = 2;
elseif(length(u) == 4)
  a(1) = u(3);
  a(2) = u(4);
  a(3) = 2;
end


f   =   @(x)([-x(1)+a(1)*(1-x(1))*exp(x(2));
		       -x(2)+a(1)*a(2)*(1-x(1))*exp(x(2))-a(3)*x(2)]);
           

% Plot the vector field, so we can directly (try) find trivial equilibrium poits : [0 1 0] and [0 0 1]
%

ax = x(1)-1;
bx = x(1)+1;

ay = x(2)-1;
by = x(2)+1;

[X,Y] = meshgrid(ax:0.1:bx , ay:0.1:by);
for i=1:1:size(X,1) 
  for j=1:1:size(X,2) 
    abu = f([X(i,j) ; Y(i,j)]);
    px(i,j) = abu(1);
    py(i,j) = abu(2);
    
    if( norm([px(i,j) py(i,j)]) ~=0)
      Npx(i,j) = px(i,j)/norm([px(i,j) py(i,j)]);
      Npy(i,j) = py(i,j)/norm([px(i,j) py(i,j)]);
    end
    
      
  end
end


figure(1);
quiver(X,Y,px,py);

quiver(X,Y,Npx,Npy);

figure;
quiver(X,Y,px,py);

hold on;
plot(x(1),x(2),'ko');

while(1)

  % wait for input
  
  [x,y,BUTTON] = ginput(1);
  
  if(BUTTON == 3)
    break;
  end

  
  
  if(BUTTON == 1)

    x0 = [ x ; y]
  
    % solve in time
    
    odefun = @(t,x)(f(x));
   
    [t,Y] = ODE45(odefun,[0:0.01:2],x0) 

    % plot iterates 
    hold on;
    plot(Y(:,1),Y(:,2),'r-');
    
  end
  if(BUTTON == 2)

    x0 = [ Y(length(Y),1) ; Y(length(Y),2) ]
  
    % solve in time
    
    odefun = @(t,x)(f(x));
   
    [t,Y] = ODE45(odefun,[0:0.01:2],x0) 

    % plot iterates 
    hold on;
    plot(Y(:,1),Y(:,2),'r-');

  end
  
  
end

