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


function tht1a(c)


% Define function f and it's derivatives

f  =  @(x,c)([ (x(2)-x(2)^2-x(1)*( x(1)^2-x(2)^2+2*x(2)^3/3 + c));
	       (x(1)+(x(2)-x(2)^2)*(x(1)^2-x(2)^2+2*x(2)^3/3+c))]);

Dxf = @(x,c)([ -(x(1)^2-x(2)^2+2*x(2)^3/3+c)-2*x(1)^2, 1-2*x(2) - x(1)*(-2*x(2)+2*x(2)^2);
	      1 + (x(2) - x(2)^2)*2*x(1), (1 - 2*x(2))*(x(1)^2-x(2)^2+2*x(2)^3/3+c)+(x(2)-x(2)^2)*(-2*x(2)+2*x(2)^2)]);
	      
Dcf = @(x,c)( [ -x(1) ; x(2)-x(2)^2]);


% Define derivative respect to u

fu  = @(u)( f(u(1:2),u(3)));
Dfu = @(u)( [Dxf(u(1:2),u(3)) Dcf(u(1:2),u(3))] );


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

[X,Y] = meshgrid(-1:0.1:1.5 , -1:0.1:1.5);
for i=1:1:size(X,1) 
  for j=1:1:size(X,2) 
    abu = f([X(i,j) ; Y(i,j)],c);
    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);


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,c));
   
    [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,c));
   
    [t,Y] = ODE45(odefun,[0:0.01:2],x0) 

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

  end
  
  
end

