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;

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


I never get a chaotic plot.

Source of the code: GitHub code

  • $\begingroup$ are you certain that the initial condition is one that leads to chaotic behavior $\endgroup$ Feb 18, 2020 at 16:47
  • $\begingroup$ What is the differential equation (or system of equations) you are trying to solve with this code? What if you replaced these equations by simple harmonic motion? Would the code work as expected? Did you look over the code and do you understand what it is doing? $\endgroup$
    – mjw
    Feb 18, 2020 at 16:58
  • $\begingroup$ @phdmba7of12 it doesn't depend only on initial condition but also on other costants (thah honestly i don't found on this code and i don't know how to implement), but i only found this code online to apply rk to van der pol in matlab $\endgroup$
    – rickamc
    Feb 18, 2020 at 16:58
  • $\begingroup$ Please find a book, such as "Numerical Recipes" by Press, et al. and look up fourth-order Runge-Kutta methods. $\endgroup$
    – mjw
    Feb 18, 2020 at 17:02
  • $\begingroup$ Looks like the code first finds $y_{1,i+1}$ from $y_{1,i}$ and then finds $y_{2,i+1}$ from both $y_{1,i}$ and $y_{2,i}$. Is that how it is supposed to work? $\endgroup$
    – mjw
    Feb 18, 2020 at 17:04

1 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;


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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .