ode45octave.m

function [T, Y] = ode45octave(FUN, tspan, Y0, tol, varargin)
%ode45octave
%   [T, Y] = ode45octave(FUN, tspan, Y0)
%
%   Integra un sistema de ecuaciones diferenciales ordinarias usando el
%   m�todo de Runge-Kutta en la variante Dormand-Price.
%
%   [T, Y] = ode45octave(FUN, tspan, Y0, tol, par_1, par_2, ...)
%      FUN      funci�n con las derivadas
%      tspan    intervalo de integraci�n [tinicial, tfinal]
%      Y0       magnitud inicial de las variables
%      tol      margen de error (por defecto 1e-10)
%      par_i    par�metros para FUN
%      T        puntos de integraci�n
%      Y        soluci�n del sistema de ecuaciones en cada punto
%
%   Esta funci�n fue creada originalmente por Marc Compere
%   http://users.powernet.co.uk/kienzle/octave/matcompat/scripts/ode_v1.11/ode45.m
%   y modificada para simplificarla y garantizar la compatibilidad con
%   Matlab. FUN debe usar la sintaxis de odefile en Matlab, de modo
%   que dYdt = fun(t, Y, [], par_1, par_2, ...)

if nargin < 4, tol = []; end
if isempty(tol), tol = 1e-10; end

pow = 1/6; % see p.91 in the Ascher & Petzold reference for more infomation.

% Use the Dormand-Prince 4(5) coefficients:
a_(1,1)=0;
a_(2,1)=1/5;
a_(3,1)=3/40;       a_(3,2)=9/40;
a_(4,1)=44/45;      a_(4,2)=-56/15;      a_(4,3)=32/9;
a_(5,1)=19372/6561; a_(5,2)=-25360/2187; a_(5,3)=64448/6561; a_(5,4)=-212/729;
a_(6,1)=9017/3168;  a_(6,2)=-355/33;     a_(6,3)=46732/5247; a_(6,4)=49/176;   a_(6,5)=-5103/18656;
a_(7,1)=35/384;     a_(7,2)=0;           a_(7,3)=500/1113;   a_(7,4)=125/192;  a_(7,5)=-2187/6784;  a_(7,6)=11/84;
% 4th order b-coefficients
b4_(1,1)=5179/57600; b4_(2,1)=0; b4_(3,1)=7571/16695; b4_(4,1)=393/640; b4_(5,1)=-92097/339200; b4_(6,1)=187/2100; b4_(7,1)=1/40;
% 5th order b-coefficients
b5_(1,1)=35/384;     b5_(2,1)=0; b5_(3,1)=500/1113;   b5_(4,1)=125/192; b5_(5,1)=-2187/6784;    b5_(6,1)=11/84;    b5_(7,1)=0;

for i=1:7
    c_(i)=sum(a_(i,:));
end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t) / 2.5;
hmin = (tfinal - t) / 1e9;
h = (tfinal - t)/100;       % initial guess at a step size
x = Y0(:);                  % this always creates a column vector, x
T = t;                      % first output time
Y = x';                     % first output solution
k_ = zeros(size(x,1), 7);   % k_ needs to be initialized as an Nx7 matrix where N=number of rows in x
% (just for speed so octave doesn't need to allocate more memory at each stage value)

% Compute the first stage prior to the main loop.  This is part of the Dormand-Prince pair caveat.
% Normally, during the main loop the last stage for x(k) is the first stage for computing x(k+1).
% So, the very first integration step requires 7 function evaluations, then each subsequent step
% 6 function evaluations because the first stage is simply assigned from the last step's last stage.
% note: you can see this by the last element in c_ is 1.0, thus t+c_(7)*h = t+h, ergo, the next step.

k_(:,1) = feval(FUN, t, x, [], varargin{1:nargin-4}); % first stage

% The main loop using Dormand-Prince 4(5) pair
while (t < tfinal) * (h >= hmin),
    if t + h > tfinal, h = tfinal - t; end
    
    % Compute the slopes by computing the k(:,j+1)'th column based on the previous k(:,1:j) columns
    % notes: k_ needs to end up as an Nxs, a_ is 7x6, which is s by (s-1),
    %        s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages)
    for j = 1:6
        k_(:,j+1) = feval(FUN, t+c_(j+1)*h, x+h*k_(:,1:j)*a_(j+1,1:j)', [], varargin{1:nargin-4});
    end
    
    % compute the 4th order estimate
    x4 = x + h * (k_*b4_);     % k_ is Nxs (or Nx7) and b4_ is a 7x1
    
    % compute the 5th order estimate
    x5 = x + h * (k_*b5_);     % k_ is Nxs (or Nx7) and b5_ is a 7x1
    
    % estimate the local truncation error
    gamma1 = x5 - x4;
    
    % Estimate the error and the acceptable error
    delta = norm(gamma1, 'inf');             % actual error
    tau = tol * max(norm(x, 'inf'), 1.0);    % allowable error
    
    % Update the solution only if the error is acceptable
    if (delta <= tau)
        t = t + h;
        x = x5;        % <-- using the higher order estimate is called 'local extrapolation'
        T = [T; t];
        Y = [Y; x'];
    end
    
    % Update the step size
    if (delta == 0.0), delta = 1e-16; end
    h = min(hmax, 0.8*h*(tau/delta)^pow);
    
    % Assign the last stage for x(k) as the first stage for computing x(k+1).
    % This is part of the Dormand-Prince pair caveat.
    % k_(:,7) has already been computed, so use it instead of recomputing it
    % again as k_(:,1) during the next step.
    k_(:,1) = k_(:,7);
    
end

if (t < tfinal), disp('Step size grew too small.'); end

if nargout == 0
    plot(T, Y, 'o-');
    T = [];
end