[SOLVED] Comparison of ODE solvers

Hi,

I was reading about Runge-Kutta methods and wanted to do a quick comparison of Euler's method, fourth order Runge-Kutta and MATLAB's ODE45.

I used the following differential equation:

yprime = @(t,y) vanderpoldemo(t,y,1);

Click here for more info on vanderpoldemo.

I was expecting the Runge-Kutta to be closer to ODE45 than Euler, but for some reason it is not. As can be seen from the plots, the Runge-Kutta solution is further away form ODE45 (which I'm assuming is the most correct one) than Euler's.

By changing the constant term in Runge-Kutta from (1/6) to about (1/3.8) I get quite close to ODE45. Here's a copy-paste of Runge-Kutta from MATLAB so you know which constant I'm talking about;

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

http://i163.photobucket.com/albums/t...ison_sys_1.jpg

http://i163.photobucket.com/albums/t...ison_sys_2.jpg

Should I just assume that my implementation of Runge-Kutta is wrong? I can't find any errors.. Here's the implementation:

Code:

`function [t,y] = runge_kutta_4order(yprime, tspan, y0, h)`

t0 = tspan(1); tfinal = tspan(end);

% set up the t values at which we will approximate the solution

t = [t0:h:tfinal]';

% include tfinal even if h does not evenly divide tfinal-t0

if t(end)~=tfinal, t=[t tfinal]; end

m = length(t);

y = [y0 zeros(length(y0), length(t)-1)];

for i=1:(m-1)

k1 = feval(yprime,t(i),y(:,i));

k2 = feval(yprime,t(i)+0.5*h, y(:,i)+0.5*h.*k1);

k3 = feval(yprime,t(i)+0.5*h, y(:,i)*0.5*h.*k2);

k4 = feval(yprime,t(i)+h, y(:,i)+h.*k3);

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

end

All comments are welcome, thanks!