%% Van der Pol Oscillator
% Zipeng, 2018-01-16
clc; clear all; close all;
% Command Window
prompt = 'Enter 1 for evaluating the effect of initial conditons, 2 for effect of \epsilon: ';
Input = input(prompt);
% Define the parameter
epsilon = [0.1 1 2];
% initial condition set
ic = {[5;0],[1;0],[-5;-5]};
% Time interval
tspan = [0:0.1:50];
switch Input
%
case 1
% Define the function
f = @(t,x) [x(2);-x(1)-epsilon(1)*(-1+x(1)^2)*x(2)];
% number of ics
[~,n] = size(ic);
% color set
col = {'r','k','b'};
%
for i = 1:n
%
[ts,ys] = ode45(f,tspan,ic{i});
hold on
plot(ys(:,1),ys(:,2),col{i},'LineWidth',1.5,'DisplayName',['x0 = ' mat2str(ic{i})]);
legend('-DynamicLegend'); % hold the legend from previous plot
%
end
%
grid on;
hold off
axis tight;
xlabel('x1');
ylabel('x2');
title('Trajectory vs Initial Conditions (\epsilon = 0.1)');
%
case 2
% number of epsilons
[~,n] = size(epsilon);
% color set
col = {'r','k','b'};
%
for i = 1:n
% Define the function
f = @(t,x) [x(2);-x(1)-epsilon(i)*(-1+x(1)^2)*x(2)];
[ts,ys] = ode45(f,tspan,ic{1});
hold on
plot(ys(:,1),ys(:,2),col{i},'LineWidth',1.5,'DisplayName',['\epsilon = ' num2str(epsilon(i))]);
legend('-DynamicLegend'); % hold the legend from previous plot
%
end
%
grid on;
hold off
axis tight;
xlabel('x1');
ylabel('x2');
title(['Trajectory vs \epsilon (x0 = ' mat2str(ic{1}) ')']);
end
%