Stability Analysis (Dertermining the Limit Cycle)

I am wanting to use a simple stability analysis technique to determine the limit cycle stability for mu< 0, for my Van Der Pol oscillator. I was ideally attempting to use a phase portrait but I don't know how to produce this in code.

The matlab code I currently have produced for the Van Der Pol oscillation can be seen below.

Any guidance would be greatly appreciated.

% Simulation parameters
DT = 0.01;                   % Time step
N  = 10000;                 % Number of discrete time points

% Equation parameter
mu = 0.1;

% Declare array to store discrete time samples
x = zeros(1,N);

% Set boundary conditions
x(1) = -0.01;
x(2) = 0.0;

% Simulate using recurrance relation
for i = 3:N   
    phi = mu*DT/2*(x(i-1)^2-1);
    x(i) = x(i-1)*(2-DT^2)/(1+phi) - x(i-2)*(1-phi)/(1+phi);
end

% Plot
plot((1:N)*DT,x,'r')

% Calculating angular frequency (omega) using period
[pks,pktimes] = findpeaks(x, (1:N)*DT);
Period = mean(diff(pktimes))
Omega = (2*pi)/Period

Answer 1

[Star Strider]

The phase portrait is usually plotted as the function against its derivative. Use the gradient function to calculate the derivative:

figure
plot(x, gradient(x))
grid

Another way is to plot the vector against a delayed version of itself:

Ofst = 5;
figure
plot(x(1:end-Ofst), x(Ofst+1:end))
grid

Both of these produce a reough phase portrait. The gradient version is likely more accurate.

[Macaulay Wright]

Thank you again! You're a real help. It worked perfectly and it was exactly what I was trying to produce.

[Star Strider]

As always, my pleasure!