%Averaging_VanderPol_Phase.m

loop=1;

xmin=-4;
xmax=4;
ymin=-4;
ymax=4;

epsilon=0.1;

%options = odeset('OutputFcn','odephas2');

figure(1)
clf;

Tfinal=input('Enter final time:  ');

while (loop==1)
   clc
   keepgoing=input('Enter 0 to end:  ');
   if keepgoing==0
      break;
   end
   ax_choice=input('Enter 0 to change axes:  ');
   if ax_choice==0
      xmin=input('xmin:  ');
      xmax=input('xmax:  ');
      ymin=input('ymin:  ');
      ymax=input('ymax:  ');
   end

   TSPAN=[0 Tfinal];
   x10=input('Enter x10:  ');
   x20=input('Enter x20:  ');
   Y0=[x10; x20];
   figure(1)
   hold on
   %[T,Y]=ode45('phase_limit_fun1',TSPAN,Y0,options);
   %[T,Y]=ode45('Averaging_VanderPol_Phase_fun',TSPAN,Y0);
   [T,Y]=ode45(@(t,y) Averaging_VanderPol_Phase_fun(t,y,epsilon),TSPAN,Y0);

   plot(Y(:,1),Y(:,2),'r-')
   axis([xmin xmax ymin ymax])
   xlabel('x1')
   ylabel('x2')
end
hold off

