## 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} =} crk(@var{dxdt}, @var{t}, @var{x0})
##
## Return the solution of differential equation @var{dxdt(t, x)} by classical Runge Kutta 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 = crk(dxdt, t, x0)
  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);
    k1 = dxdt(tn, xn);
    k2 = dxdt(tn + 0.5 * h, xn + 0.5 * h * k1);
    k3 = dxdt(tn + 0.5 * h, xn + 0.5 * h * k2);
    k4 = dxdt(tn + h, xn + h * k3);
    xn += h * to_row_vector(k1 + 2.0 * (k2 + k3) + k4) / 6.0;
    x(idx + 1, :) = xn;
  endfor
endfunction

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