% ### VDPode45EX.m ###       09.21.14
% Numerically integrate the van der Pol oscillator
%   m*x''- P.mu*(1-x^2)*x' + k*x = A*sin(wt)

clear
% -----------------------------------------------------
% User input (Note: All paramters are stored in a structure)
P.y0(1) = 0.01;   % initial position [m]
P.y0(2) = 0.0;   % initial velocity [m/s]
P.mu= 0.5;  % damping coefficient [kg/s]
P.k= 1.0;   % stiffness [N m]
P.m= 1.0;    % mass [kg]

% sinusoidal driving term
P.A= 0.0;   % amplitude [N] (set to zero to turn off)
fD= 1.05*sqrt(P.k/P.m)/(2*pi);  % freq. (Hz) [expressed as fraction of resonant freq.]

% Integration limits
P.t0 = 0.0;   % Start value
P.tf = 100.0;   % Finish value
P.dt = 0.01;  % time step
% ----------------------------------------------------------------------
% +++
% spit back out some basic derived quantities 
P.wr= 2*pi*fD;  % convert to angular freq.

% +++
% use built-in ode45 to solve
[t y] = ode45('VDPfunction', [P.t0:P.dt:P.tf],P.y0,[],P);

% ------------------------------------------------------
% visualize
figure(1); clf;
plot(t,y(:,1)); hold on; grid on; 
xlabel('t [s]');    ylabel('x(t) [m]')
% Phase plane 
figure(2); clf;
plot(y(:,1), y(:,2)); hold on; grid on; 
xlabel('x [m]');    ylabel('dx/dt [m/s]')