how can i draw Vector/Direction field for van der pol oscillator ?

I am having problem in drawing the vector / direction field for Van der Pol Oscillator which looks something like this https://en.wikipedia.org/wiki/Van_der_Pol_oscillator#/media/File:VanDerPolOscillator.png

i have written this code :

options = odeset('MaxStep',0.5);
temp = inputdlg('Enter mu value');
mu = str2double(temp{1,1});

[t,y] = ode45(@(t,y) vdp1_1(t,y,mu),[0 10],[2; 0],options);

figure(1);
plot(t,y(:,1)), 
xlabel('time t');
ylabel('Y1(t)');
title('t vs Y1(t)');

figure(2);
plot(t,y(:,2))
xlabel('time t');
ylabel('Y2(t)');
title('t vs Y2(t)');

figure(3),
 plot(y(:,1),y(:,2)) , hold on
quiver(y(:,1), y(:,2), gradient(y(:,1)), gradient(y(:,2) ))
hold off
xlabel('y1(t)');
 ylabel('Y2(t)');
 title('y1(t) vs Y2(t)');

function dydt = vdp1_1(t,y,mu)

dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = [mu * (1-y(1)^2)*y(2)-y(1)];

end

my code doesnot give any errors , and also doesnt give the vector field .It just gives gradient boundary of the ode using quiver for y1 vs y2

Any help would be appreciable!

[KSSV]

Same question: https://in.mathworks.com/matlabcentral/answers/331449-how-to-draw-the-direction-field-of-van-der-pol-oscillator

Answer 1

[KSSV]

You have to follow like shown in the link:

http://matlab.cheme.cmu.edu/2011/08/09/phase-portraits-of-a-system-of-odes/

Answer 2

[Michael Mathew]

function VanderPol_plot()
 
steps = 20;
 
s1 = linspace(-6,6,steps);
s2 = linspace(-6,6,steps);
 
% creates two matrices one for all the x-values on the grid, and one for
% all the y-values on the grid.
 
[x,y] = meshgrid(s1,s2);
 
u = zeros(size(x));
v = zeros(size(x));
 
for i=1:numel(x)
    [t,y_sol] = ode23(@vdp1,[0 20],[x(i); y(i)]);
    figure(1);
    plot(y_sol(:,1),y_sol(:,2)); hold on;
    figure(2);
    quiver(y_sol(:,1), y_sol(:,2), gradient(y_sol(:,1)), gradient(y_sol(:,2) ), 'r'); hold on;
end
 
figure(1);
ylabel('Solution plot');
figure(2);
title('Gradient plot')
 
end
function dydt = vdp1(t,y, mu)
%Vandepol oscillator ode
if nargin < 3
    mu = 1;
end
dydt = [y(2); (mu-y(1)^2)*y(2)-y(1)];
end