% clean up left-behind garbage
clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Draw a phase plane of the Van der Pol oscillator. %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This program can be modified to draw phase planes of two-dimensional
% first order autonomous systems, of the form
%    y1' = f1(y1,y2)
%    y2' = f2(y1,y2)
% where primes indicate differentiation to some time-like parameter t.

% Simple scalar first order ordinary differential equations of the form
% d(y)/d(x) = f(x,y) can also be handled by setting y1 =x, y2 = y, and f1
% and f2 so that their ratio equals f.

% To modify the program, you will need to create a function file, call it
% myfunc.m, that returns f1 and f2.  Pattern myfunc.m after file
% vanderpol.m. Then in the code below, search for vanderpol and replace it
% with myfunc.  You may also need to make changes in the program setup
% variables that follow now.

% Program setup variables:

% plot limits
y1min = -5;
y1max = 5;
y2min = -5;
y2max = 5;

% number of direction-field arrows in each direction
y1cnt = 20;
y2cnt = 20;

% arrow scaling factor
arrscale = .35;

% maximum Runge-Kutta step count
stepmax = 1000;

% To test for accuracy, reduce the next variables by a factor 2 and see
% whether it makes a difference:

% maximum Runge-Kutta relative step size
dyrmax = .03;

% maximum Runge-Kutta angular change
dangmax = .1;

% maximum Runge-Kutta magnitude change factor
dmagmax = 2;

% minimum Runge-Kutta relative step size before we call it a day
dyrmin = .0001;

% End of the program setup variables.

% compact output please
format compact

% create the direction field
y1pos = linspace(y1min,y1max,y1cnt+2);
y2pos = linspace(y2min,y2max,y2cnt+2);
[y1m,y2m] = meshgrid(y1pos(2:y1cnt+1),y2pos(2:y2cnt+1));
f1m = zeros(size(y1m));
f2m = zeros(size(y1m));
for i = 1:y2cnt
  for j = 1:y1cnt
    k1 = vanderpol([y1m(i,j);y2m(i,j)]);
    k11 = k1(1);
    k12 = k1(2);
    k1m = sqrt(k11^2+k12^2);
    if k1m == 0
      k1m = 1;
    end
    f1m(i,j) = k11/k1m;
    f2m(i,j) = k12/k1m;
  end
end
quiver(y1m,y2m,f1m,f2m,arrscale,'r')
xlabel('y_1')
ylabel('y_2')
axis([y1min y1max y2min y2max], "equal")
hold on
fig = gcf()
set (fig, 'numbertitle', 'off', 'name', 'Phase Plane: Click a point to plot or press q to quit.')

% loop over the mouse clicks
while 1 == 1

% get the initial conditions from the user's mouse click
[y10, y20, button] = ginput(1)
if button > 1
  break
end

% solve the ODE in both directions
for sgn = -1:2:1

  % initialize
  t = 0;
  y = [y10; y20];

  % first plot point
  clear y1 y2
  i = 1;
  y1(i) = y(1);
  y2(i) = y(2);
  if y(1) < y1min || y(1) > y1max || y(2) < y2min || y(2) > y2max
    break
  end

  % initialize the first time step
  k1 = sgn*vanderpol(y);
  k11 = k1(1)/((y1max-y1min));
  k12 = k1(2)/((y2max-y2min));
  k1m = sqrt(k11^2+k12^2);
  invdt =  max(abs(k11),abs(k12));
  if invdt == 0
    break
  end
  dt = dyrmax/invdt;

  for step = 1:stepmax

    % initialize the time step, allowing 50 % growth
    k1 = sgn*vanderpol(y);
    k11 = k1(1)/((y1max-y1min));
    k12 = k1(2)/((y2max-y2min));
    k1m = sqrt(k11^2+k12^2);
    invdt =  max(abs(k11),abs(k12));
    if invdt == 0
      break
    end
    dt0 = dt;
    dt = dyrmax/invdt;
    if dt > 1.5*dt0
      dt = 1.5*dt0;
    end

    % reduce the step size as needed
    ok = 0;
    dt = dt*2;
    while ok ~= 1
      dt = dt/2;
      ok = 1;
      k2 = sgn*vanderpol(y + .5*k1*dt);
      k21 = k2(1)/((y1max-y1min));
      k22 = k2(2)/((y2max-y2min));
      k2m = sqrt(k21^2+k22^2);
      if k2m > k1m*(1+.5*(dmagmax-1))
	ok = 0;
      end
      if k2m < k1m/(1+.5*(dmagmax-1))
	ok = 0;
      end
      if k11*k21+k12*k22 < 0
	ok = 0;
      end
      if abs(k11*k22-k12*k21) > dangmax*k1m*k2m
	ok = 0;
      end
      if dt == 0
	ok = 1;
      end
    end

    % do the time step
    k2 = sgn*vanderpol(y + .5*k1*dt);
    k3 = sgn*vanderpol(y + .5*k2*dt);
    k4 = sgn*vanderpol(y + k3*dt);
    y = y + (k1 + 2*k2 + 2*k3 + k4)*dt/6;

    % store the plot point
    i = i+1;
    y1(i) = y(1);
    y2(i) = y(2);

    if y(1) < y1min || y(1) > y1max || y(2) < y2min || y(2) > y2max
      break
    end

    if k1m*dt < dyrmin
      break
    end

  end

  i
  if i > 1
    plot(y1,y2)
  end

end

end
