Hi there,

ODE45 is, essentially, an RK4 method, and so you should see almost no error between the results that ode45 give, and the results that RK4 give.

I just ran my own code, with my own rk4 method, and the from rk4 and ode45 are identical, almost no difference whatsoever.

With that said, I can't see any particular problem with your implementation.

So, what I'll do is give you my code, and allow you to spot any errors yourself.

It could be that your stepsize is too large? What stepsize are you using?

Here is my code. I haven't commented it, but hopefully it is self-explanatory.

Code:

clear
clc
StartTime = 0;
EndTime = 20;
Stepsize = 0.01;
x = [2,2];
xdot = [0,0];
tout=zeros(1,(EndTime/Stepsize)+1);
xout=zeros((EndTime/Stepsize)+1, length(x));
xdout = xout;
i=0;
for time = StartTime:Stepsize:EndTime,
i = i+1 ;
tout(i) = time ;
xout(i,:) = x ;
xdout(i,:) = xdot ;
xdot = vanderpoldemo(time,x,1);
x = rk4('vanderpoldemo',Stepsize,time,x);
end
[TOUT,YOUT] = ode45(@vanderpoldemo,[StartTime EndTime],[2;2]);
figure(1)
hold on
plot(tout,xout(:,1),'r--')
plot(TOUT,YOUT(:,1),'b--')
title('Comparison of ODE Solvers')
xlabel('time (s)')
ylabel('Function value')
legend('Runge-Kutta 4','ODE45')
figure(2)
hold on
plot(tout,xout(:,2),'r--')
plot(TOUT,YOUT(:,2),'b--')
title('Comparison of ODE Solvers')
xlabel('time (s)')
ylabel('Function value')
legend('Runge-Kutta 4','ODE45')

Code:

function x = rk4(model,h,t,x)
k1 = feval(model,t,x,1) ;
k2 = feval(model,t+0.5*h,x+h*k1'/2,1) ;
k3 = feval(model,t+0.5*h,x+h*k2'/2,1) ;
k4 = feval(model,t+h,x+h*k3',1) ;
x = x + h*(k1' + 2*k2' + 2*k3' + k4')/6 ;
end

Code:

function dydt = vanderpoldemo(t,y,Mu)
Mu = 1;
dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

Let me know how you get on.