%=======================================================================
% Chapter 9, code 2: adaptive RKF45 for n first-order ODEs.
% Based on the revised Fortran companion code.
% Mathematical example:
%   x1' = x2, x2' = -x1, x3' = -2*x3, x4' = x1+x3.
% Initial values at t = 0: (1,0,1,0).
% Exact: (cos(t), -sin(t), exp(-2*t), sin(t)+(1-exp(-2*t))/2).
%
% To supply another system:
%   1. Set n and all initial values in the driver.
%   2. Edit equations: dx(j) is the derivative of x(j).
%   3. Choose the interval, first trial step, and tolerance.
%   4. Disable check_exact or edit the exact_solution function.
% All adaptive steps are inside rkf45_system; it needs no changes.
% State vectors are columns. Returned X has one time per row and one
% component per column: X(i,j) is x_j at t(i).
%=======================================================================
clear;

%-------------------- Initial data: edit here ----------------------------
n = 4;
t0 = 0;
t_end = 10;
x0 = zeros(n,1);
x0(1) = 1;
x0(2) = 0;
x0(3) = 1;
x0(4) = 0;
h_first = 0.5;
tolerance = 1.0e-8;      % Same absolute local tolerance per component
check_exact = true;
%------------------------------------------------------------------------
[t,X,h,accepted,rejected] = ...
    rkf45_system(@equations,n,t0,x0,t_end,h_first,tolerance);

% Print and save the initial point and every accepted point.
filename = 'ch09_rkf45_system.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open data file: %s',filename);
end
fprintf('Adaptive RKF45: system of first-order ODEs\n');
fprintf('Number of equations: %d\n',n);
fprintf('Tolerance: %.4e\n',tolerance);
fprintf('%10s','t');
fprintf(fid,'# t');
for j = 1:n
    fprintf(' %18s',sprintf('x%d',j));
    fprintf(fid,' x%d',j);
end
fprintf(' %18s\n','h');
fprintf(fid,' h\n');
for i = 1:length(t)
    fprintf('%10.6f',t(i));
    fprintf(' %18.10e',X(i,:));
    fprintf(' %18.10e\n',h(i));
    fprintf(fid,'%.16e',t(i));
    fprintf(fid,' %.16e',X(i,:));
    fprintf(fid,' %.16e\n',h(i));
end
fclose(fid);
fprintf('Accepted steps: %8d\n',accepted);
fprintf('Rejected trials:%8d\n',rejected);
if check_exact
    % This exact solution matches ONLY the default system and initial data.
    exact = exact_solution(t_end);
    fprintf('%10s %18s %18s %12s\n', ...
            'component','numerical x','exact x','abs error');
    for j = 1:n
        fprintf('%10d %18.10e %18.10e %12.4e\n', ...
                j,X(end,j),exact(j),abs(X(end,j)-exact(j)));
    end
end
fprintf('Data file: %s\n',filename);

% One panel per component; panel count follows n.
% Exact curves use a dense grid; markers show actual accepted points.
if check_exact
    t_plot = linspace(t0,t_end,1000)';
    X_plot = exact_solution(t_plot);
end
figure;
for j = 1:n
    subplot(n,1,j);
    if check_exact
        plot(t_plot,X_plot(:,j),'k-','LineWidth',1.2);
        hold on;
        plot(t,X(:,j),'ko--','MarkerSize',3);
        hold off;
        if j == 1
            legend('Exact','RKF45','Location','best');
        end
    else
        plot(t,X(:,j),'ko-','MarkerSize',3);
    end
    ylabel(sprintf('x_%d(t)',j));
    grid on;
    if j == 1
        title('Adaptive RKF45: system of first-order ODEs');
    end
    if j == n
        xlabel('t');
    end
end

% h(i) is the accepted step used to reach t(i); omit the initial h = 0.
figure;
semilogy(t(2:end),h(2:end),'ko-','MarkerSize',3);
xlabel('t');
ylabel('Accepted step size h');
title('Adaptive RKF45: step size');
grid on;

%==================== Equations and exact solution ======================
function dx = equations(t,x)
% One derivative per line, in the same order as the initial values.
% Return a column vector of ALL derivatives; do not change x.
% t is available for equations that depend explicitly on time.
dx = zeros(size(x));
dx(1) = x(2);
dx(2) = -x(1);
dx(3) = -2*x(3);
dx(4) = x(1)+x(3);
end

function X = exact_solution(t)
% For the DEFAULT system only. One time per row, one component per column.
t = t(:);
X = [cos(t), -sin(t), exp(-2*t), sin(t)+(1-exp(-2*t))/2];
end

%==================== Numerical method ==================================
function [t_values,X,h_values,accepted,rejected] = ...
    rkf45_system(f,n,t0,x0,t_end,h_first,tolerance)
% Adaptive Fehlberg RKF45; advance with the fifth-order estimate.
% The largest absolute embedded component error controls acceptance.
% One tolerance assumes comparable numerical scales; rescale variables
% or use component tolerances for widely differing magnitudes or units.
% Local error control does not bound accumulated global error.
max_attempts = 100000;
h_min = 1.0e-12;
if n < 1 || length(x0) ~= n || t_end <= t0 || h_first <= 0 || tolerance <= 0
    error('Use n initial values, t_end > t0, and positive step and tolerance.');
end
t = t0;
x = x0(:);
h = h_first;
accepted = 0;
rejected = 0;
attempts = 0;
t_values = t0;
X = x';
h_values = 0;

while t < t_end
    if t_end-t <= 64*eps(1)*max(1,abs(t_end))
        break;
    end
    if attempts >= max_attempts
        error('Maximum number of trials exceeded.');
    end
    h = min(h,t_end-t);
    tf = t+h;
    if h < h_min || tf <= t
        error('Step too small to advance time.');
    end
    attempts = attempts+1;

    % Each k is a column vector and includes h.
    % All equations use the SAME complete intermediate state at each stage.
    k1 = h*f(t,x);
    k2 = h*f(t+h/4,x+k1/4);
    k3 = h*f(t+3*h/8,x+3*k1/32+9*k2/32);
    k4 = h*f(t+12*h/13,x+1932*k1/2197-7200*k2/2197+7296*k3/2197);
    k5 = h*f(t+h,x+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
    k6 = h*f(t+h/2,x-8*k1/27+2*k2-3544*k3/2565 ...
               +1859*k4/4104-11*k5/40);

    trial = x+16*k1/135+6656*k3/12825+28561*k4/56430 ...
             -9*k5/50+2*k6/55;
    component_error = abs(k1/360-128*k3/4275-2197*k4/75240+k5/50+2*k6/55);
    err = max(component_error);
    if err <= tolerance
        if h == t_end-t
            tf = t_end;
        end
        t = tf;
        x = trial;
        accepted = accepted+1;
        t_values(end+1,1) = t;
        X(end+1,:) = x';
        h_values(end+1,1) = h;
    else
        rejected = rejected+1; % Keep both t and x unchanged.
    end

    % Embedded error scales approximately as h^5.
    if err == 0
        factor = 2;
    else
        factor = 0.9*(tolerance/err)^0.2;
    end
    if err <= tolerance
        h = h*min(2,factor);
    else
        h = h*max(0.2,factor);
    end
end
end