function [tout, yout, ireject] = stdrk75(fcn, gcn, t0, tfinal, y0, tol)
% SDRK75: Second-derivative explicit Runge-Kutta method (7th/5th order)
% Solves y' = f(y), y'' = f'(y) * f(y) = g(y)
% Inputs:
%   fcn  - function handle for f(y)
%   gcn  - function handle for g(y)
%   t0   - initial time
%   tfinal - final time
%   y0   - initial condition
%   tol  - tolerance for adaptive step size
% Outputs:
%   tout - time values
%   yout - solution values
%   ireject - number of rejected steps

% Coefficients
a = [  0          0         0        0       0       0;
      1/98        0         0        0       0       0;
     -1/98     5/49         0        0       0       0;
    169/1024 -119/2048 357/2048      0       0       0;
    -29/18   231/85   -112/135   512/2295    0       0;
     11/270  2401/12240 2401/12960 512/6885 1/288    0]';

b  = [11/270, 2401/12240, 2401/12960, 512/6885, 1/288, 0];
bb = [53/270, -(343/2448), 6517/12960, -832/6885, -11/288, 1/10];
c  = [0, 1/7, 3/7, 3/4, 1, 1]';

% Initialization
t = t0;
y = y0(:);
s = length(c);
g = zeros(length(y), s);

% Estimate initial step size
f = fcn(y);
pow = 1/7;
h = tol^pow / max(max(abs(f)), 1e-2);
hmax = (tfinal - t) / 5;
hmin = (tfinal - t) / 2e6;
h = min(hmax, max(h, hmin));

% Preallocate output buffers (initial guess)
max_steps = ceil((tfinal - t0) / hmin);
tout = zeros(max_steps, 1);
yout = zeros(max_steps, length(y));
tout(1) = t;
yout(1, :) = y.';
step = 1;

% Initial second derivative
g(:,1) = gcn(y);
ireject = 0;

% Main integration loop
while t < tfinal && h >= hmin
    if t + h > tfinal
        h = tfinal - t;
    end

    % Compute g at intermediate stages
    for j = 2:s
        y_stage = y + c(j)*h*f + h^2 * g * a(:,j);
        g(:,j) = gcn(y_stage);
    end

    % Error estimation (embedded method)
    delta = norm(h * g * (b - bb)', inf)^(1.1666);

    if delta <= tol
        % Accept step
        t = t + h;
        y = y + h * f + h^2 * g * b';
        f = fcn(y);               % next step's f
        g(:,1) = g(:,s);          % FSAL: reuse last g
        step = step + 1;
        tout(step) = t;
        yout(step, :) = y.';
    else
        ireject = ireject + 1;
    end

    % Adjust step size
    if delta ~= 0
        h = min(hmax, 0.8 * h * (tol / delta)^pow);
    end
end

% Trim unused preallocated space
tout = tout(1:step);
yout = yout(1:step, :);

% Warn if failed
if t < tfinal
    warning('Singularity likely. Integration stopped at t = %g.', t);
end
end


%---------------------------------------------------------------------------
% A sample run for Kaps problem with $\xi=200$ and $\delta=10^{-9}$ could be

>> xi=200;
[tout, yout, ireject] = sdrk75(@(y) [-y(1)*(1+y(1))+y(2); ...
                                xi*(y(1)^2-y(2))-2*y(2)], ...
    @(y) [y(1)+(3 + xi)*y(1)^2+2*y(1)^3-(xi + 3)*y(2)-2*y(1)*y(2);    ...
    -(4*xi+xi^2)*y(1)^2-2*xi*y(1)^3+2*xi*y(1)*y(2)+(xi + 2)^2*y(2)],  ...
    0, 10*pi, [1 1]',  1e-9);
fprintf('stages=%6.0f\n', length(tout)*6 + ireject*5);
fprintf('max abs error=%9.2e\n',max(max(abs(yout - [exp(-tout) exp(-2*tout)]))));

stages= 11073
max abs error= 7.72e-10