Code snippet: Runge-Kutta-Gil (Matlab)

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;