%%%%%%%%%
% run.m %
%%%%%%%%%
h = 0.02;
MIN_X = -6;
MAX_X =  2;
MIN_Y = -3;
MAX_Y =  3;

clf;
more off;
axis ([MIN_X MAX_X MIN_Y MAX_Y]);
grid on;
hold on;

button = 0;
while (1)
 if button==0  [x0 y0 button] = ginput(1); end;
 if (button==3) break; end;

 plot(x0,y0,'r.');
 over_bound_forward = 0;
 over_bound_backward = 0;
 x_last_forward = x0;
 y_last_forward = y0;
 x_last_backward = x0;
 y_last_backward = y0;

 while(1)

  if ~over_bound_forward 
  x = [x_last_forward y_last_forward];
  %%%%%%%%%%%%
  % Forward  %
  %%%%%%%%%%%%
  for i=1:1000,
   xt(i) = x(1);
   yt(i) = x(2);
   % 4-th order Runge-Kutta integration.
   k1 = F(x)*h;
   k2 = F(x+k1/2.)*h;
   k3 = F(x+k2/2.)*h;
   k4 = F(x+k3)*h;
   x = x+(k1+k2+k2+k3+k3+k4)/6.;
   if (x(1)<MIN_X)|(x(1)>MAX_X)|(x(2)<MIN_Y)|(x(2)>MAX_Y) 
    over_bound_forward = 1; 
    fprintf (1,'Forward overbound: xlast = %f ylast = %f \n',x(1),x(2));
    break; 
   end;
  end;
  x_last_forward = xt(i);
  y_last_forward = yt(i);
  plot(xt(1:i),yt(1:i),'--');
  axis ([MIN_X MAX_X MIN_Y MAX_Y]);
  end;


  if ~over_bound_backward 
  x = [x_last_backward y_last_backward];
  %%%%%%%%%%%%
  % Backward %
  %%%%%%%%%%%%
  for i=1:1000,
   xt(i) = x(1);
   yt(i) = x(2);
   % 4-th order Runge-Kutta integration.
   k1 = F(x)*h;
   k2 = F(x-k1/2.)*h;
   k3 = F(x-k2/2.)*h;
   k4 = F(x-k3)*h;
   x = x-(k1+k2+k2+k3+k3+k4)/6.;
   if (x(1)<MIN_X)|(x(1)>MAX_X)|(x(2)<MIN_Y)|(x(2)>MAX_Y) 
    over_bound_backward = 1; 
    fprintf (1,'Backward overbound: xlast = %f ylast = %f \n',x(1),x(2));
    break; 
   end;
  end;
  x_last_backward = xt(i);
  y_last_backward = yt(i);
  plot(xt(1:i),yt(1:i));
  axis ([MIN_X MAX_X MIN_Y MAX_Y]);
  end;

  if over_bound_forward & over_bound_backward 
  over_bound_forward = 0;
  over_bound_backward = 0;
  button = 0;
  break; 
  end;

  fprintf(1, 'Continue from previous trajectory? If YES press middle button.. \n');
  [x1 y1 button] = ginput(1);
  if (button~=2) x0 = x1;  y0 = y1; fprintf(1,'This trail finished. \n\n'); break; end;
 end;
end;

xlabel('X');
ylabel('Y');
axis ([MIN_X MAX_X MIN_Y MAX_Y]);
ss = input('The phase portrait of ','s');
title(['Phase Portrait of ' ss]);

PRINT;
more on; 



%%%%%%%%%
% run.m %
%%%%%%%%%
function [f] = F(X)
x = X(1);
y = X(2);
y2 = y*y;
f(1) = (-x-(y2-1))*y2*y;
f(2) = (y2-1)*((2+x)+(3-y2));
end;