function [] = VanderPol()
% VanderPol.m
% MST 383/683 
% Fall 2020
%
% Function calculating and plotting the solution of the Van der Pol 
% oscillator using ode45
%
% Code compiled from the ode45 documentation:
% https://www.mathworks.com/help/matlab/ref/ode45.html

% Test out ode45 using a single ODE:
tspan = [0 5]; % initial time = 0, final time = 5
y0 = 0;        % initial condition
[t,yy] = ode45(@(t,y) 2*t*y, tspan, y0); % dy/dt = 2t

plot(t,yy,'-o')

% Now approximate the solution to Van der Pol using ode45:
clear

mu = 1;

tspan = [0 20];
z0 = [2; 0];    % two ICs needed for a 2nd order ODE [x0, y0]
[t,z] = ode45(@vdp1,tspan,z0);

figure, hold on
set(gca,'FontSize',12)
plot(t,z(:,1),'-o','LineWidth',2)
plot(t,z(:,2),'-o','LineWidth',2)
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating a function inside another function creates a ~nested~ function. 
%
% Variable values are automatically passed to nested functions - they turn
% blue to indicate the value is shared across multiple functions.

function f = vdp1(t,z)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

x = z(1);
y = z(2);

f = [y; mu*(1-x^2)*y-x];
end

end