User:MikeSkyHigh/'''Matlab'''/Nonlinear Control Course Code/Assignment 1/Q5 Van del Pol Oscillator

%% 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      
%