RK4 for Van der Pol MatLab

Question

I'm trying to implement RK algorithm applied to Van Der Pol equation. I tried this code

y = zeros(n, 2);
y(1, 1) = y0(1);
y(1, 2) = y0(2);

% Van der Pol equation 
% y''(t) - (1 - y^2(t)) * y'(t) + y(t) = 0;
% y'(t)  = y(2)
% y''(t) = (1-y(1)^2)*y(2)-y(1)

for i = 1:n-1
    k1 = h*(y(i,2));
    k2 = h*(y(i,2)+0.5*k1);
    k3 = h*(y(i,2)+0.5*k2);
    k4 = h*(y(i,2)+k3);
    y(i+1,1) = y(i,1)+(k1+2*k2+2*k3+k4)/6;

    y1 = y(i,1);
    y2 = y(i,2);
    k1 = h*((1-y1^2)*y2-y1);

    y1 = y(i,1)+0.5*k1;
    y2 = y(i,2)+0.5*k1;
    k2 = h*((1-y1^2)*y2-y1);

    y1 = y(i,1)+0.5*k2;
    y2 = y(i,2)+0.5*k2;
    k3 = h*((1-y1^2)*y2-y1);

    y1 = y(i,1)+k3;
    y2 = y(i,2)+k3;
    k4 = h*((1-y1^2)*y2-y1);

    y(i+1,2) = y(i,2)+(k1+2*k2+2*k3+k4)/6;
end

Trying to give different h, t_end, y0 as input but with

plot(y)

I never get a chaotic plot.

Source of the code: GitHub code

Answer 1

You need to implement RK4 code for the the coupled system as a coupled system. You can avoid the code bloat that would usually imply by utilizing the vector capabilities of the MATrix LABoratory

But staying with the component-wise code you should get

y1 = y(i,1); v1 = y(i,2)
k1y = h*v1;
k1v = h*((1-y1^2)*v1 - y1);

y2 = y1 + 0.5*k1y; v2 = v1 + 0.5*k1v;
k2y = h*v2;
k2v = h*((1-y2^2)*v2 - y2);

y3 = y1 + 0.5*k2y; v3 = v1 + 0.5*k2v;
k3y = h*v3;
k3v = h*((1-y3^2)*v3 - y3);

y4 = y1 + k3y; v4 = v1 + k3v;
k4y = h*v4;
k4v = h*((1-y4^2)*v4 - y4);

y(i+1,1) = y(i,1)+(k1y+2*k2y+2*k3y+k4y)/6;
y(i+1,2) = y(i,2)+(k1v+2*k2v+2*k3v+k4v)/6;

etc.

Note that the Van der Pol oscillator is an oscillator, that it, is has a limit cycle which is a deformation of the circle of radius $2$. If you want a halfway chaotic plot, you should implement the Lorenz system.