## Copyright (C) 2010, 2016, 2020  Shinji Tamura  <stamura@ecs.shimane-u.ac.jp>
##
## This program 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 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{x} =} rkf(@var{dxdt}, @var{t}, @var{x0})
##
## Return the solution of differential equation @var{dxdt(t, x)} by Runge Kutta Fehlberg formula.
##
## @var{x} is the matrix, with each row and column corresponding to the elements
## of vector @var{t} and state vector, respectively.
##
## Arguments:
## @itemize
##   @item @var{dxdt} : function representing the differential equation,
##       1st argument is time(scalar) and 2nd argument is state vector
##   @item @var{t} : time vector
##   @item @var{x0} : initial value vector for @var{x}
## @end itemize
##
## @end deftypefn

function x = rkf(dxdt, t, x0)
  C = [0.0, 0.25, 3.0 / 8.0, 12.0 / 13.0, 1.0, 0.5];
  A1 = 1.0 / 4.0;
  A2 = [3.0, 9.0] / 32.0;
  A3 = [1932.0, -7200.0, 7296.0] / 2197.0;
  A4 = [8341.0, -32832.0, 29440.0, -845.0] / 4104.0;
  A5 = [-6080.0, 41040.0, -28352.0, 9295.0, -5643.0] / 20520.0;
  B = [902880, 0.0, 3953664.0, 3855735.0, -1371249.0, 277020.0] / 7618050.0;
  ##
  if (rows(x0) > columns(x0))
    to_row_vector = @(x) transpose(x);
  else
    to_row_vector = @(x) x;
  endif
  x = zeros(numel(t), numel(x0));
  x(1, :) = to_row_vector(x0);
  xn = to_row_vector(x0);
  x1 = dxdt(t(1), x0);
  if (rows(x1) > columns(x1))
    to_row_vector = @(x) transpose(x);
  else
    to_row_vector = @(x) x;
  endif
  ##
  for idx = 1 : numel(t) - 1
    h = t(idx + 1) - t(idx);
    tn = t(idx) + h * C;
    k1 = to_row_vector(dxdt(tn(1), xn));
    k2 = to_row_vector(dxdt(tn(2), xn + h * A1 * k1));
    k3 = to_row_vector(dxdt(tn(3), xn + h * A2 * [k1; k2]));
    k4 = to_row_vector(dxdt(tn(4), xn + h * A3 * [k1; k2; k3]));
    k5 = to_row_vector(dxdt(tn(5), xn + h * A4 * [k1; k2; k3; k4]));
    k6 = to_row_vector(dxdt(tn(6), xn + h * A5 * [k1; k2; k3; k4; k5]));
    xn += h * B * [k1; k2; k3; k4; k5; k6];
    x(idx + 1, :) = xn;
  endfor
endfunction

##
## Local Variables:
## mode: octave
## End:
