%
%  vdp.m -- 07 nov 03 -- djm
%
%  - van der pol oscillator
%  - ralf's friday problem

clear
global ODE ep mu

%  set ODE choice
disp(' ')
disp(' ODE choices:')
disp('   1) van der pol')
disp('   2) vdp (lienard)')
disp(' ')
setup = input('choose #/R:  ');
disp(' ')
disp(' click to choose IV (top to continue; bottom to quit)')
disp(' ')

%  initialization
switch setup
    case{1}
        ODE = 1;
        ep = 1;  mu = 4.0;  box = 8.0;
        pname = '\bf van der pol oscillator';
        yname = '\bf x''(t)-axis';
    case{2}
        ODE = 2;
        ep = 1;  mu = 10.0;  box = 4.0;
        pname = '\bf vdp -- lienard form';
        yname = '\bf y(t)-axis';
    otherwise
        break
end
   
%  initial condition
figure(1);  clf;  hold on
if (ODE==2)
    xp = -box:box/100:box;
    plot(xp,(xp.^3)/3 - xp,'b')
end
title(pname)
xlabel('\bf x(t)-axis')
ylabel('\bf x''(t)-axis')
axis(box*[-1 1 -1 1]);  
text(box/4, 1.1*box,'\bf click to continue -> X')
text(box/4,-1.15*box,'\bf click to quit -> X')

IV = [];

%  plot & solve loop
while (isempty(IV))
    click = ginput(1);
    IV = round(10*click)/10;
   
%  bottom to quit, top to continue
    
    if(IV(2) < -box)
        disp(['IV = quit'])
        break
    elseif(IV(2) > box)
        IV = y(end,:);
        disp(['IV = continue'])
    else
        disp(['IV = ' num2str(IV)])
    end
    
%  invoke runge-kutta solver
    
    opt = odeset('RelTol',1e-3,'AbsTol',1e-6);
    [t y] = ode45(@Fvdp,[0 50],IV,opt);

    plot(y(:,1),y(:,2),'r')
    plot(IV(1),IV(2),'ko')
    
    IV = [];
end

figure(2);  clf
[t y] = ode45(@Fvdp,[0 100],IV,opt);
plot(t,y(:,1),'b')
hold on
plot(t,y(:,2),'b--')
title(['\bf final IV plots:  \mu = ' num2str(mu)])
xlabel('\bf t-axis')
ylabel('\bf x(t) = solid & x''(t) or y(t) = dashed')