From CFD-Wiki
function [tout,xout] = myrkg(fname,tspan,x0)
%
% Based on code written by:
% Marc Compere
% CompereM@asme.org
% created : 06 October 1999
% modified: 17 January 2001
%
% Copyright (C) 2001, 2000 Marc Compere
%
% rkg.m is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% <krivilli at unberwoot.net> - Modified to use Runge-Kutta-Gil method (06/03/2006):
% without adaptive step, possible stiffness related issues.
%
% Kunge-Kutta-Gil coefficients
w(1)=1/6;
w(2)=(2-sqrt(2))/6;
w(3)=(2+sqrt(2))/6;
w(4)=1/6;
c(1)=0;
c(2)=1/2;
c(3)=1/2;
c(4)=1;
a(2,1)=1/2;
a(3,1)=-(sqrt(2)-1)/2;
a(3,2)=(2-sqrt(2))/2;
a(4,1)=0;
a(4,2)=-sqrt(2)/2;
a(4,3)=(2+sqrt(2))/2;
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
h = (tfinal - t)/500.0; % '500' is our step. Change it as needed.
x = x0(:); % column vector
tout = t;
xout = x.';
while (t < tfinal)
if t + h > tfinal, h = tfinal - t; end
ydot=feval(fname,t,x);
k(:,1)=h*ydot;
ydot=feval(fname,t+c(2)*h,x+a(2,1)*k(:,1));
k(:,2)=h*ydot;
ydot=feval(fname,t+c(3)*h,x+a(3,1)*k(:,1)+a(3,2)*k(:,2));
k(:,3)=h*ydot;
ydot=feval(fname,t+c(4)*h,x+a(4,2)*k(:,2)+a(4,3)*k(:,3));
k(:,4)=h*ydot;
x2=x+w(1)*k(:,1)+w(2)*k(:,2)+w(3)*k(:,3)+w(4)*k(:,4);
t = t + h;
x = x2;
tout = [tout; t];
xout = [xout; x.'];
end;