function [tout, yout, dyout, varargout] = rkn1210(funfcn, tspan, y0, yp0, options)
% RKN1210       12th/10th order Runge-Kutta-Nystrom integrator
% 
% RKN1210() is a 12th/10th order numerical integrator for ordinary
% differential equations of the special form
% 
%   y'' = f(t, y)                     (1)
% 
% with initial conditions
% 
%   y (t0) = y0                       (2)
%   y'(t0) = yp0 
% 
% This second-order differential equation is integrated with a
% Runge-Kutta-Nystrom method, with 17 function evaluations per step. The
% RKN-class of integrators is especially suited for this purpose, since
% compared to a classic Runge-Kutta integration scheme the same accuracy
% can be obtained with less function evaluations. 
% 
% This RKN12(10) method is a very high-order method, to be used in problems 
% with *extremely* stringent error tolerances. As the name implies, the
% error should be less than O(h^13). In verious studies, it has been shown 
% that this particular integration technique is overall more efficient for 
% ODE's of the form (1) than multi-step or extrapolation methods that give 
% the same accuracy.
% 
% RKN1210()'s behavior is very similar MATLAB's ODE-integrator suite:
% 
% USAGE: 
% ----------------
% 
%   [t, y, yp] = RKN1210(funfcn, tspan, y0, yp0)
%   [t, y, yp] = RKN1210(funfcn, tspan, y0, yp0, options)
%
%   [t, y, yp, exitflag, output] = RKN1210(...)
%   [t, y, yp, TE, YE, YPE, IE, exitflag, output] = RKN1210(...)
% 
% INPUT ARGUMENTS:
% ----------------
% 
%  funfcn  - definition of the second-derivative function f(t, y) 
%            (See (1)). It should accept a scalar value [t] and a 
%            column vector [y] which has the same number of elements 
%            as the initial values [y0] and [dy0] provided. 
% 
%    tspan - time interval over which to integrate. It can be a
%            two-element vector, in which case it will be interpreted
%            as an interval. In case [tspan] has more than two 
%            elements, the integration is carried out to all times in
%            [tspan]. Only the values for those times are then
%            returned in [y] and [yp].
% 
%  y0, yp0 - initial values as in (2). Both should contain the same number
%            of elements.
% 
%  options - options structure, created with ODESET(). Used options are
%            MaxStep, InitialStep, AbsTol, Stats, Event, OutputFcn, 
%            OutputSel, and Refine. See the help for ODESET() for more 
%            information.
%
% How to use Event and/or Output functions is described in the documentation
% on ODESET(). There is one difference: RKN1210() now also passes the first 
% derivative [yp] at each step as an argument:
%
%   status = outputFcn(t, y, yp, flag)
%   [value, isterminal, direction] = event(t, y, yp)
%
% where [t] is scalar, and [y] and [yp] are column vectors, as with f(t,y).
% 
%
% OUTPUT ARGUMENTS:
% ----------------
% 
% t, y, yp - The approximate solutions for [y] and [y'] at times [t].
%            All are concatenated row-wise, that is 
%
%               t  = N-by-1
%               y  = N-by-numel(y0)
%               y' = N-by-numel(y0)
%
%            with N the number of sucessful steps taken during the 
%            integration. 
% 
%  exitflag - A scalar value, indicating the termination conditions 
%             of the integration:
%                 
%             -2: a non-finite function value was encountered during the
%                 integration (INF of NaN); the integration was stopped.
%             -1: the step size [h] fell below  the minimum acceptable 
%                 value at some time(s) [t]; results may be inaccurate.
%              0: nothing was done; initial state.
%             +1: sucessful integration, normal exit.
%             +2: integration was stopped by one of the output
%                 functions.
%             +3: One or more events were detected, and their 
%                 corresponding [isterminal] condition also evaluated to 
%                 [true]. 
%
% TE,YE,    - These arguments are only returned when one or more event
%  YPE,IE     functions are used. [TE] contains the times at which events
%             were detected. [YE] and [YPE] lists the corresponding values 
%             of the solution [y] and the first derivative [yp] at these 
%             times. [IE] contains indices to the event-functions with 
%             which these events were detected. Use a smaller value for
%             AbsTol (in [options]) to increase the accuracy of these 
%             roots when required.
% 
%    output - structure containing additional information about the
%             integration. It has the fields:
%
%             output.h              step size (sucesful steps only) at 
%                                   each time [tn] 
%             output.rejected       amount of rejected steps
%             output.accepted       amount of accepted steps
%             output.delta          estimate of the largest possible 
%                                   error at each time [tn]   
%             output.message        Short message describing the
%                                   termination conditions
% 
%             Note that these fields contain the information of ALL
%             steps taken, even for cases where [tspan] contains
%             more than 2 elements. 
% 
%
% See also ODE45, ODE86, RKN86.



% Author:
% Name       : Rody P.S. Oldenhuis
% E-mail     : oldenhuis@gmail.com
% Affiliation: LuxSpace sarl.
%
% please report any bugs or suggestions to oldenhuis@gmail.com


    
% Based on the code for ODE86 and RKN86, also available on the MATLAB
% FileExchange. 
% 
% The construction of RKN12(10) is described in  
% High-Order Embedded Runge-Kutta-Nystrom Formulae
% J. R. DORMAND, M. E. A. EL-MIKKAWY, AND P. J. PRINCE
% IMA Journal of Numerical Analysis (1987) 7, 423-430
% 
% Coefficients obtained from 
% http://www.tampa.phys.ucl.ac.uk/rmat/test/rknint.f
% These are also available in any format on request to these authors.  



%{
ELEMENTARY EXAMPLE 
(2-body gravitational interaction / circular orbit):

clc

% Equations of motion for a circular orbit in 2D
f1 = @(t, y) -y/sqrt(y'*y)^3;
f2 = @(t, y) [y(3:4); -y(1:2)/sqrt(y(1:2).'*y(1:2))^3];

% RKN1210 with moderate accuracy setting  
disp('RKN1210, with AbsTol = RelTol = 1e-6:')
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
tic
  [t1, y1] = rkn1210(f1, [0 1000], [1; 0], [0; 1], options);
toc  
disp('Maximum absolute error:')
max( (y1(:,1)-cos(t1)).^2 + (y1(:,2)-sin(t1)).^2 )

% This is how much ODE45 will have to be tuned to achieve similar accuracy
fprintf('\n\n')  
disp('RKN1210, with AbsTol = RelTol = 1e-6:')
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
tic
  [t2, y2] = ode45(f2, [0, 1000], [1; 0; 0; 1], options);
toc  
disp('Maximum absolute error:')
max( (y2(:,1)-cos(t2)).^2 + (y2(:,2)-sin(t2)).^2 )

%}
  
  
%{
CHANGELOG: 
    2013/June/11 
       - Set 'InitialStep' to previous step on recursive calls to RKN1210

    2013/June/07
       - Added value of the second derivative for succesful steps to the output structure
       - No memory is assigned for output variables when no output arguments are requested       
       - Removed all support for funtions defined be string; only function handles 
       - Added a few more checks on user inputs
       - Updated elementary example
       - Minor stylistic changes
%}

    %% Error traps
    
    % args error check
    narginchk(4,5); 
    
    % further checks on input 
    assert( isa(funfcn, 'function_handle'),...
        'rkn1210:funfcn_isnt_a_function', ...
        'Second derivative f(t,y) must be given as function handle.');
    assert( ~isempty(tspan) && numel(tspan) >= 2,...
        'rkn1210:tspan_empty', 'Time interval [tspan] must contain at least 2 values.');
    assert( all(diff(tspan) ~= 0),...
        'rkn1210:tspan_dont_span', 'Values in [tspan] must all span a non-zero interval.');
    assert( all(diff(tspan) > 0) || all(diff(tspan) < 0), ...
        'rkn1210:tspan_must_increase', 'The entries in tspan must strictly increase or decrease.');    
    assert( numel(y0) == numel(yp0),...
        'rkn1210:initial_values_disagree', ...
        'Initial values [y0] and [yp0] must contain the same number of elements.');
    if nargin == 5
        assert( isstruct(options),...
            'rkn1210:options_not_struct', ...
            'Options must be given as a structure created with ODESET().');
    end
    
    %% Load the coefficients
    
    % c
    c = [0.0e0
        2.0e-2
        4.0e-2
        1.0e-1
        1.33333333333333333333333333333e-1
        1.6e-1
        5.0e-2
        2.0e-1
        2.5e-1
        3.33333333333333333333333333333e-1
        5.0e-1
        5.55555555555555555555555555556e-1
        7.5e-1
        8.57142857142857142857142857143e-1
        9.45216222272014340129957427739e-1
        1.0e0
        1.0e0];
    
    % matrix A is lower triangular. It's easiest to 
    % load the coefficients row-by-row:
    A = zeros(17);
    A(2,1:1) = 2.0e-4;  
    
    A(3,1:2) = [2.66666666666666666666666666667e-4   
                5.33333333333333333333333333333e-4];    
            
    A(4,1:3) = [2.91666666666666666666666666667e-3
                -4.16666666666666666666666666667e-3
                6.25e-3];    
            
    A(5,1:4) = [1.64609053497942386831275720165e-3
                0.0e0
                5.48696844993141289437585733882e-3
                1.75582990397805212620027434842e-3];    
            
    A(6,1:5) = [1.9456e-3
                0.0e0
                7.15174603174603174603174603175e-3
                2.91271111111111111111111111111e-3
                7.89942857142857142857142857143e-4];    
            
    A(7,1:6) = [5.6640625e-4
                0.0e0
                8.80973048941798941798941798942e-4
                -4.36921296296296296296296296296e-4
                3.39006696428571428571428571429e-4
                -9.94646990740740740740740740741e-5];    
            
    A(8,1:7) = [3.08333333333333333333333333333e-3
                0.0e0
                0.0e0
                1.77777777777777777777777777778e-3
                2.7e-3
                1.57828282828282828282828282828e-3
                1.08606060606060606060606060606e-2];    
            
    A(9,1:8) = [3.65183937480112971375119150338e-3
                0.0e0
                3.96517171407234306617557289807e-3
                3.19725826293062822350093426091e-3
                8.22146730685543536968701883401e-3
                -1.31309269595723798362013884863e-3
                9.77158696806486781562609494147e-3
                3.75576906923283379487932641079e-3];
    
    A(10,1:9) = [3.70724106871850081019565530521e-3
                0.0e0
                5.08204585455528598076108163479e-3
                1.17470800217541204473569104943e-3
                -2.11476299151269914996229766362e-2
                6.01046369810788081222573525136e-2
                2.01057347685061881846748708777e-2
                -2.83507501229335808430366774368e-2
                1.48795689185819327555905582479e-2];
    
    A(11,1:10) = [3.51253765607334415311308293052e-2
                0.0e0
                -8.61574919513847910340576078545e-3
                -5.79144805100791652167632252471e-3
                1.94555482378261584239438810411e0
                -3.43512386745651359636787167574e0
                -1.09307011074752217583892572001e-1
                2.3496383118995166394320161088e0
                -7.56009408687022978027190729778e-1
                1.09528972221569264246502018618e-1];
    
    A(12,1:11) = [2.05277925374824966509720571672e-2
                0.0e0
                -7.28644676448017991778247943149e-3
                -2.11535560796184024069259562549e-3
                9.27580796872352224256768033235e-1
                -1.65228248442573667907302673325e0
                -2.10795630056865698191914366913e-2
                1.20653643262078715447708832536e0
                -4.13714477001066141324662463645e-1
                9.07987398280965375956795739516e-2
                5.35555260053398504916870658215e-3];
    
    A(13,1:12) = [-1.43240788755455150458921091632e-1
                0.0e0
                1.25287037730918172778464480231e-2
                6.82601916396982712868112411737e-3
                -4.79955539557438726550216254291e0
                5.69862504395194143379169794156e0
                7.55343036952364522249444028716e-1
                -1.27554878582810837175400796542e-1
                -1.96059260511173843289133255423e0
                9.18560905663526240976234285341e-1
                -2.38800855052844310534827013402e-1
                1.59110813572342155138740170963e-1];
    
    A(14,1:13) = [8.04501920552048948697230778134e-1
                0.0e0
                -1.66585270670112451778516268261e-2
                -2.1415834042629734811731437191e-2
                1.68272359289624658702009353564e1
                -1.11728353571760979267882984241e1
                -3.37715929722632374148856475521e0
                -1.52433266553608456461817682939e1
                1.71798357382154165620247684026e1
                -5.43771923982399464535413738556e0
                1.38786716183646557551256778839e0
                -5.92582773265281165347677029181e-1
                2.96038731712973527961592794552e-2];
    
    A(15,1:14) = [-9.13296766697358082096250482648e-1
                0.0e0
                2.41127257578051783924489946102e-3
                1.76581226938617419820698839226e-2
                -1.48516497797203838246128557088e1
                2.15897086700457560030782161561e0
                3.99791558311787990115282754337e0
                2.84341518002322318984542514988e1
                -2.52593643549415984378843352235e1
                7.7338785423622373655340014114e0
                -1.8913028948478674610382580129e0
                1.00148450702247178036685959248e0
                4.64119959910905190510518247052e-3
                1.12187550221489570339750499063e-2];
    
    A(16,1:15) = [-2.75196297205593938206065227039e-1
                0.0e0
                3.66118887791549201342293285553e-2
                9.7895196882315626246509967162e-3
                -1.2293062345886210304214726509e1
                1.42072264539379026942929665966e1
                1.58664769067895368322481964272e0
                2.45777353275959454390324346975e0
                -8.93519369440327190552259086374e0
                4.37367273161340694839327077512e0
                -1.83471817654494916304344410264e0
                1.15920852890614912078083198373e0
                -1.72902531653839221518003422953e-2
                1.93259779044607666727649875324e-2
                5.20444293755499311184926401526e-3];
    
    A(17,1:16) = [1.30763918474040575879994562983e0
                0.0e0
                1.73641091897458418670879991296e-2
                -1.8544456454265795024362115588e-2
                1.48115220328677268968478356223e1
                9.38317630848247090787922177126e0
                -5.2284261999445422541474024553e0
                -4.89512805258476508040093482743e1
                3.82970960343379225625583875836e1
                -1.05873813369759797091619037505e1
                2.43323043762262763585119618787e0
                -1.04534060425754442848652456513e0
                7.17732095086725945198184857508e-2
                2.16221097080827826905505320027e-3
                7.00959575960251423699282781988e-3
                0.0e0];   
            
    % this facilitates and speeds up implementation       
    A = A.';            
    
    % Bhat (high-order b)
    Bhat = [1.21278685171854149768890395495e-2
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            8.62974625156887444363792274411e-2
            2.52546958118714719432343449316e-1
            -1.97418679932682303358307954886e-1
            2.03186919078972590809261561009e-1
            -2.07758080777149166121933554691e-2
            1.09678048745020136250111237823e-1
            3.80651325264665057344878719105e-2
            1.16340688043242296440927709215e-2
            4.65802970402487868693615238455e-3
            0.0e0
            0.0e0];
    
    % BprimeHat (high-order b-prime)
    Bphat  = [1.21278685171854149768890395495e-2
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            9.08394342270407836172412920433e-2
            3.15683697648393399290429311645e-1
            -2.63224906576909737811077273181e-1
            3.04780378618458886213892341513e-1
            -4.15516161554298332243867109382e-2
            2.46775609676295306562750285101e-1
            1.52260530105866022937951487642e-1
            8.14384816302696075086493964505e-2
            8.50257119389081128008018326881e-2
            -9.15518963007796287314100251351e-3
            2.5e-2];
    
    % B (low-order b)
    B = [1.70087019070069917527544646189e-2
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            7.22593359308314069488600038463e-2
            3.72026177326753045388210502067e-1
            -4.01821145009303521439340233863e-1
            3.35455068301351666696584034896e-1
            -1.31306501075331808430281840783e-1
            1.89431906616048652722659836455e-1
            2.68408020400290479053691655806e-2
            1.63056656059179238935180933102e-2
            3.79998835669659456166597387323e-3
            0.0e0
            0.0e0];
    
    % Bprime (low-order bprime)
    Bp = [1.70087019070069917527544646189e-2
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            0.0e0
            7.60624588745593757356421093119e-2
            4.65032721658441306735263127583e-1
            -5.35761526679071361919120311817e-1
            5.03182602452027500044876052344e-1
            -2.62613002150663616860563681567e-1
            4.26221789886109468625984632024e-1
            1.07363208160116191621476662322e-1
            1.14139659241425467254626653171e-1
            6.93633866500486770090602920091e-2
            2.0e-2
            0.0e0];
    
    %% Initialize
    
    % initialize all variables
       exitflag = 0;                          pow = 1/12;    
             t0 = tspan(1);                tfinal = tspan(end);
              t = t0;                           y = y0(:);
             dy = yp0(:);                    tout = t0;
           yout = y0(:).';                  dyout = yp0(:).';    
           hmin = abs(tfinal - t)/1e12;        f  = y0(:)*zeros(1,17);
    have_events = false;                     halt = false;
 have_outputFcn = false;           produce_output = (nargout ~= 0);
    
    % initialize output-structure
           output.h = [];               output.fevals = 0;
    output.rejected = 0;              output.accepted = 0;
       output.delta = [];             output.message  = 'Integration not started.';
      output.d2ydt2 = [];
                  
    % times might be given in reverse
    direction = 1 - 2*(numel(tspan)==2 && tfinal<t0);       
   
    % parse options
    if (nargin < 5)
        options = odeset; end                             % initial options structure
    
    Stats       = odeget(options, 'Stats', 'off');        % display statistics at the end?
    abstol      = odeget(options, 'AbsTol', 1e-14);       % Absolute tolerance
    reltol      = odeget(options, 'RelTol', 1e-7);        % Relative tolerance
    hmax        = odeget(options, 'MaxStep', tfinal-t);   % defaults to entire interval
    initialstep = odeget(options, 'InitialStep');         % default determined later
    Event       = odeget(options, 'Event', []);           % defaults to no eventfunctions
    OutputFcn   = odeget(options, 'OutputFcn', []);       % defaults to no output functions    
    
    
    % FUTURE WORK    
    NormControl = odeget(options, 'NormControl', 'off');  % relax error control
    MassMatrix  = odeget(options, 'Mass', []);            % Mass matrix, for problems of the ...
    % (and friends)                                       % form M(t,y)*y'' = F(t,y)
    
    
    % in case of Event/output functions, define a few extra parameters    
    if ~isempty(Event)
        
        % so we DO have to evaluate event-functions
        have_events = true;
        % cast into cell-array if only one function is given (easier later on)
        if isa(Event, 'function_handle')
            Event = {Event}; end
        % number of functions provided
        num_events = numel(Event); 
        
        % Check if all are indeed function handles
        for ii = 1:num_events
            if isempty(Event{ii}) || ~isa(Event{ii}, 'function_handle')
                error('rkn1210:event_not_function_handles',...
                    ['Unsupported class for event function received; event %d is of class ''%s''.\n',...
                    'RKN1210 only supports function handles.'], ii, class(Event{ii}));
            end
        end
        
        % initialize TE (event times), YE (event solutions) YPE (event
        % derivs) and IE (indices to corresponding event function). Check
        % user-provided event functions at the same time
        previous_event_values = zeros(num_events,1);
        for k = 1:num_events
            try
                previous_event_values(k) = Event{k}(t, y, dy); 
                
            catch ME
                ME2 = MException('rkn1210:eventFcn_dont_evaluate',...
                    sprintf('Event function #%1d failed to evaluate on initial call.', k));
                rethrow(addCause(ME,ME2));  
                
            end
        end
        
        if produce_output        
            TE = []; YE = []; YPE = []; IE = []; end        
        
    end
    
    if ~isempty(OutputFcn)
        
        % so we DO have to evaluate output-functions
        have_outputFcn = true;
        % cast into cell-array if only one function is given (easier later on)
        if isa(OutputFcn, 'function_handle')
            OutputFcn = {OutputFcn}; end
        % number of functions provided
        num_outputFcn = numel(OutputFcn);     
        
        % Check if all are indeed function handles
        for ii = 1:num_outputFcn
            if isempty(OutputFcn{ii}) || ~isa(OutputFcn{ii}, 'function_handle')
                error('rkn1210:output_not_function_handles',...
                    ['Unsupported class for output function received; output function %d is of class ''%s''.\n',...
                    'RKN1210 only supports function handles.'], ii, class(OutputFcn{ii}));
            end
        end
        
        % WHICH elements should be passed to the output function?
        OutputSel = odeget(options, 'OutputSel', 1:numel(y0)); 
        % adjust the number of points passed to the output functions by this factor
        Refine = odeget(options, 'Refine', 1);         
                
        % call each output function with 'init' flag. Also check whether
        % the user-provided output function evaluates
        for k = 1:num_outputFcn
            try
                OutputFcn{k}(t, y(OutputSel), dy(OutputSel), 'init');
                
            catch ME
                ME2 = MException('rkn1210:OutputFcn_doesnt_evaluate',...
                    sprintf('Output function #%d failed to evaluate on initial call.', k));
                rethrow(addCause(ME,ME2));                
                
            end
        end
    end
    
    %% Different use case: numel(tspan) > 2
    
    % do the calculation recursively if [tspan] has 
    % more than two elements
    if (numel(tspan) > 2)
        
        % the output times are already known
        tout = tspan(:);
        
        % call this function as many times as there are times in [tspan]
        for ii = 1:numel(tspan)-1
            
            % new initial values
            tspanI = tspan(ii:ii+1);
            y0     = yout(end, :);
            yp0    = dyout(end, :);
                        
            % call the integrator
            if have_events
                [toutI, youtI, dyoutI, TEI, YEI, DYEI, IEI, exitflag, outputI] = ...
                    rkn1210(funfcn, tspanI, y0, yp0, options);                
            else
                [toutI, youtI, dyoutI, exitflag, outputI] = ...
                    rkn1210(funfcn, tspanI, y0, yp0, options);
            end
            
            % new initial step is old next-to-last step
            options = odeset(options, ...
                'InitialStep', outputI.h(end-1)); 
            
            % append the solutions
            yout  = [yout;  youtI(end, :)];  %#ok
            dyout = [dyout; dyoutI(end, :)]; %#ok
            if have_events
                TE   = [TE; TEI];  %#ok
                YEI  = [YE; YEI];  %#ok
                DYEI = [YPE; DYEI];%#ok
                IEI  = [IE; IEI];  %#ok
            end
            
            % process the output
            output.h        = [output.h; outputI.h];
            output.fevals   = output.fevals + outputI.fevals;
            output.rejected = output.rejected + outputI.rejected;
            output.accepted = output.accepted + outputI.accepted;
            output.delta    = [output.delta; outputI.delta];
            
            % evaluate any output functions at each [t] in [tspan]
            if have_outputFcn
                % evaluate all functions
                for k = 1:num_outputFcn
                    try
                        halt = OutputFcn{k}(toutI, youtI(OutputSel), dyoutI(OutputSel), []);
                        
                    catch ME
                        ME2 = MException('rkn1210:OutputFcn_failure_integration',...
                            sprintf('Output function #%1d failed to evaluate during integration.', k));
                        rethrow(addCause(ME,ME2));
                        
                    end
                end
                % halt integration when requested
                if halt
                    exitflag = 2; 
                    finalize(); 
                    return
                end
            end
            
            % should we quit?
            if exitflag == -2 || exitflag == 3
                break, end
        end
        
        % we're done.
        index = size(yout,1); 
        finalize(); 
        return;
        
    end % if, for more than 2 elements in tspan

    %% Initial step
    
    try
        f(:,1) = funfcn(t, y);
        
    catch ME
        ME2 = MException('rkn1210:incorrect_funfcnoutput', sprintf(...
            'Derivative function should return a %3.0f-element column vector.', numel(y0)));
        rethrow(addCause(ME,ME2));
        
    end
    
    output.fevals = output.fevals + 1; % don't forget this one :)
    
    if isempty(initialstep)
        % default initial step
        h = abstol^pow / max(max(abs([dy.' f(:, 1).'])), 1e-4);
        h = min(hmax,max(h,hmin));        
        
    else
        % user provided initial step
        h = initialstep;
        
    end
    
    h  = direction*h; % take care of direction
    h2 = h*h; % pre-compute the square
    
    %% The main loop
    
    % Initialize all variables
    index = 1;
    grow_arrays();
    
    % the main loop    
    while ( abs(t-tfinal) > 0)
        
        % take care of final step
        if ((t + h) > tfinal)
            h  = (tfinal - t); 
            h2 = h*h; % also pre-compute the square
        end

        % Compute the second-derivative
        % NOTE: 'Vectorized' in ODESET() has no use; we need the UPDATED
        % function values to calculate the NEW ones, i.e., the function
        % evaluations are not independent. 
        for jj = 1:17
            f(:,jj) = funfcn( t + c(jj)*h, y + c(jj)*h*dy + h2*f*A(:,jj) );            
            output.fevals = output.fevals + 1;
        end
        
        % check for inf or NaN
        if any(~isfinite(f(:)))
            exitflag = -2;
            % use warning (not error) to preserve output thus far
            warning('rkn1210:nonfinite_values',...
                 ['INF or NAN value encountered during the integration.\n',...
                  'Terminating integration...']);   
            finalize();
            return;
        end % non-finite values

        % pre-compute the sums of the products with the coefficients
        fBphat = f*Bphat;   fBhat  = f*Bhat;

        % Estimate the error and the acceptable error
        delta1 = max(abs(h2*(fBhat - f*B))); % error ~ |Y - y|
        delta2 = max(abs(h*(fBphat - f*Bp)));% error ~ |dot{Y} - dot{y}|
        delta  = max(delta1, delta2);        % worst case error

        % update the solution only if the error is acceptable
        if (delta <= abstol) &&...
           (delta/max(norm(y+h*dy+h2*fBhat),norm(dy+h*fBphat)) <= reltol)
            
            % update the new solution
            index = index + 1;
            tp    = t;
            t     = t + h;  
            yp    = y;
            y     = y + h*dy + h2*fBhat;
            dyp   = dy;
            dy    = dy + h*fBphat; 
            
            % if new values won't fit, first grow the arrays
            % NOTE: This construction is WAY better than growing the arrays on
            % EACH iteration; especially for "cheap" integrands, this
            % construction causes a lot less overhead.
            if index > size(yout,1)
                grow_arrays(); end
            
            % insert updated values
            if produce_output
                tout (index,:) = t;
                yout (index,:) = y.';
                dyout(index,:) = dy.';
                
                output.d2ydt2(index,:) = f(:,1).';
                output.h     (index,:) = h;
                output.delta (index,:) = delta;
                output.accepted = output.accepted + 1;
            end
                                    
            % evaluate event-funtions
            if have_events
                % evaluate all event-functions
                for k = 1:num_events
                    % evaluate event function, and check if any have changed sign
                    %
                    % NOTE: although not really necessary (event functions have been 
                    % checked upon initialization), use TRY-CATCH block to produce 
                    % more useful errors in case something does go wrong. 
                    terminate = false;
                    try
                        % evaluate function
                        [value, isterminal, zerodirection] = ...
                            Event{k}(t, y, dy);
                            
                        % look for sign change                              
                        if (previous_event_values(k)*value < 0)   
                            % ZERODIRECTION: 
                            %  0: detect all zeros (default
                            % +1: detect only INcreasing zeros
                            % -1: detect only DEcreasing zeros
                            if (zerodirection == 0) ||...
                               (sign(value) == sign(zerodirection))
                                % terminate?
                                terminate = terminate || isterminal;
                                % detect the precise location of the zero
                                % NOTE: try-catch is necessary to prevent things like
                                % discontinuous event-functions from resulting in
                                % unintelligible error messages
                                try
                                    detect_Events(k, tp, previous_event_values(k), t, value);
                                    
                                catch ME
                                    ME2 = MException('rkn1210:eventFcn_failure_zero',...
                                        sprintf('Failed to locate a zero for event function #%1d.', k));
                                    throwAsCaller(addCause(ME,ME2));
                                    
                                end
                            end
                        end
                        
                        % save new value
                        previous_event_values(k) = value;
                        
                    catch ME
                        ME2 = MException('rkn1210:eventFcn_failure_integration',...
                            sprintf('Event function #%1d failed to evaluate during integration.', k));
                        rethrow(addCause(ME,ME2));
                        
                    end
                    
                    % do we need to terminate?
                    if terminate
                        exitflag = 3; 
                        finalize(); 
                        return; 
                    end
                    
                end
            end % Event functions
            
            % evaluate output functions
            if 	have_outputFcn
                % evaluate all event-functions               
                for k = 1:num_outputFcn
                    % evaluate output function
                    try
                        halt = OutputFcn{k}(t, y(OutputSel), dy(OutputSel), []);
                        
                    catch ME
                        ME2 = MException('rkn1210:OutputFcn_failure_integration',...
                            sprintf('Output function #%1d failed to evaluate during integration.', k));
                        rethrow(addCause(ME,ME2));
                        
                    end
                    % halt integration when requested
                    if halt
                        exitflag = 2;
                        finalize(); 
                        return;
                    end
                    
                end                
            end % Output functions
                    
        % rejected step: just increase its counter
        else
            output.rejected = output.rejected + 1;
        end % accept or reject step

        % adjust the step size
        if (delta ~= 0)
            h  = sign(h)*min(hmax, 0.9*abs(h)*(abstol/delta)^pow);
            % Use [Refine]-option when output functions are present
            if have_outputFcn, h = h/Refine; end
            % pre-compute the square
            h2 = h*h;
            % check the new stepsize
            if (abs(h) < hmin) 
                exitflag = -1;
                % use warning to preserve results thus far
                warning('rkn1210:stepsize_too_small', ...
                    ['Failure at time t = %6.6e: \n',...
                    'Step size fell below the minimum acceptible value of %6.6d.\n',...
                    'A singularity is likely.'], t, hmin);
            end            
        end % adjust step-size
        
    end % main loop
    
    % if the algorithm ends up here, all was ok
    finalize();
    
    %% Finalize
    
    % clean up and finalize
    function finalize
        
        % first cut off redundant zeros
        if produce_output
            tout  = tout (1:index,:);
            yout  = yout (1:index,:);
            dyout = dyout(1:index,:);
            output.h     = output.h(1:index,:);
            output.delta = output.delta(1:index,:);
        end
        
        % neutral flag means all OK
        if (exitflag == 0)
            exitflag = 1; end    
        
        % final call to the output functions
        if have_outputFcn
            for kk = 1:num_outputFcn
                try
                    OutputFcn{kk}(t, y(OutputSel), dy(OutputSel), 'done');
                    
                catch ME                    
                    ME2 = MException('rkn1210:OutputFcn_failure_finalization',...
                        sprintf('Output function #%1d failed to evaluate during final call.', k));
                    rethrow(addCause(ME,ME2));
                    
                end
            end
        end     
        
        % add message to output structure
        if produce_output
            switch exitflag
                case +1
                    output.message = 'Integration completed sucessfully.';
                case +2
                    output.message = 'Integration terminated by one of the output functions.';
                case +3
                    output.message = 'Integration terminated by one of the event functions.';
                case -1
                    output.message = sprintf(['Integration terminated successfully, but the step size\n',...
                        'fell below the minimum allowable value of %6.6e.\n',...
                        'for one or more steps. Results may be inaccurate.'], hmin);
                case -2
                    output.message = sprintf(['Integration unsuccessfull; second derivative function \n',...
                        'returned a NaN or INF value in the last step.']);
            end
                        
            % handle varargout
            if have_events
                varargout{1} = TE;
                varargout{2} = YE;
                varargout{3} = YPE;
                varargout{4} = IE;
                varargout{5} = exitflag;
                varargout{6} = output;
                
            else
                varargout{1} = exitflag;
                varargout{2} = output;
                
            end
            
        end
        
        % display stats
        if strcmpi(Stats, 'on')               
            fprintf(1, '\n\n Number of successful steps     : %6d\n', output.accepted);
            fprintf(1, ' Number of rejected steps       : %6d\n', output.rejected);
            fprintf(1, ' Number of evaluations of f(t,y): %6d\n\n', output.fevals);            
            fprintf(1, '%s\n\n', output.message);
        end
        
    end % finalize the integration 
    
    %% Detect events
    
    % Use false-position method (derivative-free)
    function detect_Events(which_event, ta, fa, tb, fb)
        
        % initialize
        opts = options;
        iterations = 0;
        maxiterations = 1e4;
        y0  = yp;
        dy0 = dyp;
        tt  = ta;
        
        % prune unnessesary options, and set initial step to current step
        opts = odeset(opts,...
            'Event'      , [],...
            'OutputFcn'  , [],...
            'Stats'      , 'off',...
            'InitialStep', h);
        
        % start iteration
        while (min(abs(fa),abs(fb)) > abstol)
            
            % make step
            iterations = iterations + 1;
            ttp = tt;
            tt  = (0.5*fb*ta - fa*tb) / (0.5*fb-fa);
            
            % safety precaution
            if (ttp == tt)
                break, end
            
            % Evaluating the event-function at this new trial location is
            % somewhat complicated. We need to recursively call this
            % RKN1210-routine to get appropriate values for [y] and [dy] at
            % the new time [tt] into the event function:
            [~, Zyout, Zdyout, ~, Zoutput] = ...
                rkn1210(funfcn, [ttp, tt], y0, dy0, opts);
            
            % set new initial step to next-to-last step of previous call
            opts = odeset(opts, 'InitialStep', Zoutput.h(end-1)); 
            
            % save old values for next iteration
            y0 =  Zyout(end,:).';
            dy0 = Zdyout(end,:).';
            
            % NOW evaluate event-function with these values
            fval = Event{which_event}(tt, y0, dy0);
            
            % keep track of number of function evaluations
            if produce_output
                output.fevals = output.fevals + Zoutput.fevals; end
            
            % compute new step
            if (fb*fval>0)
                tb = tt; fb = fval;
            else
                ta = tt; fa = fval;
            end
            
            % check no. of iterations
            if (iterations > maxiterations)
                error('rkn1210:rootfinder_exceeded_max_iterations',...
                    'Root could not be located within %d iterations. Exiting...',...
                    maxiterations);
            end
            
        end % Regula-falsi loop
            
        % The zero has been found!
        if produce_output
            
            %  insert values into proper arrays        
            TE  = [TE; tt];   YPE = [YPE; dy0];                       
            YE  = [YE; y0];   IE  = [IE; which_event];  

            % if new values won't fit, first grow the arrays
            if index > size(yout,1)
                grow_arrays(); end

            % The integrand first overshoots the zero; that's how it's
            % detected. We want the zero to be in the final arrays, but we also
            % want them in chronological order. So, move the overshoot one
            % down, and insert the zero in its place:

            index = index + 1;  

            yout (index,:) = yout (index-1,:);     yout (index-1,:) = y0.';
            dyout(index,:) = dyout(index-1,:);     dyout(index-1,:) = dy0.';
            tout (index,:) = tout (index-1,:);     tout (index-1,:) = tt;

            output.h(index,:) = tt-tout(index-1,:);
        end
        
    end % find zeros of Event-functions
    
    %% Grow arrays
    
    % the "grow" function
    function grow_arrays
        if produce_output
            
            init   = 5e2;
            nans   = NaN(init,1);
            nans_y = NaN(init,numel( y0));
            
            tout  = [tout;  nans];
            yout  = [yout;  nans_y];
            dyout = [dyout; nans_y];
            
            output.delta  = [output.delta;  nans  ];
            output.h      = [output.h;      nans  ];
            output.d2ydt2 = [output.d2ydt2; nans_y];
            
        end
    end
    
end % RKN1210 integrator
