vderpol solves the second order Van der Pol equation

function vderpol
%
% vderpol solves the second order Van der Pol equation
%
% w" -k*(1-w*w)*w' + w = 0
%
% as an autonomous first order system in v (=w') and w
% over a time interval 'dt'. It prompts for starting
% coordinates 'w0' and 'v0,' and requests
% permission to complete another time step.
% The derivatives specifying the equation are given in
% vdplfn.m
%
% vderpol.m uses: vdplfn.m, ode23k.m
%

% Initialization
disp('First order system for Van der Pol equation.')
disp(' ')
disp('When you press RETURN you will see an empty plot of v vs w.')
disp('You will be asked for the model parameter k, and initial coordinates.')
disp('Press RETURN to end any pauses.')
pause
more = 'y';
dt = 0.5; trace = 0; tol = 0.05;
xp = [-5 5]; yp = [-5 5];
k = 0;
clg; hold off
plot(xp,yp,'.')
title('1st order system for Van der Pol equation')
xlabel('w')
ylabel('v')
hold on
pause

% Input parameters
k = input(' Model parameter k : ');
while more ~= 'n'
init = input('Initial coordinates (from -4 to 4), enter as [w0 v0]: ');
dt=input('Suggest a time step (try 0.5 if you have no better suggestion): ');
disp('please wait...')
y0 = init'; t0 = 0; tf = dt; step = 'y'; mti=1;
while step == 'y'
%
[tout,yout] = ode23k('vdplfn', k, t0, tf, y0, tol, trace);
%
xp = yout(:,1); yp=yout(:,2);
nti = length(tout); mtf=mti+nti-1;
xpp(mti:mtf)=tout(1:nti);
ypp(mti:mtf)=xp(1:nti); mti=mtf+1;
plot(xp,yp,'r');
pause
disp(' t w(t) v(t)')
disp([tf xp(nti) yp(nti) ])
step = input('Another time step? (y/n): ','s');
t0 = tf; tf = t0 + dt;
num = size(tout); y0 = yout(num(1),:)';
end
more = input('Do you want another solution? (y/n): ','s');
end
hold off
clg
disp('Press ENTER for a plot of the last solution: w(t) vs t ');
pause
plot(xpp,ypp,'g');
title('Van der Pol equation')
xlabel('t')
ylabel('w')