%=======================================================================
% Chapter 9, code 1: adaptive RKF45 for one first-order ODE x' = f(t,x).
% Based on the revised Fortran companion code.
% The driver supplies initial data. All stages and adaptive step decisions
% are inside rkf45. The initial point and every accepted point are returned.
%
% Example: x' = 100/(1+10000*t^2), x(-1) = atan(-100).
% Exact: x(t) = atan(100*t). Rapid change near zero makes adaptation visible.
% To solve another equation, change the initial data and function f.
% Disable check_exact or update BOTH exact formulas (comparison and plot).
%=======================================================================
clear;

%-------------------- Initial data: edit here ----------------------------
t0 = -1;
x0 = atan(100*t0);
t_end = 1;
h_first = 0.1;
tolerance = 1.0e-6;      % Absolute local error tolerance
check_exact = true;
%------------------------------------------------------------------------
[t,x,h,accepted,rejected] = rkf45(@f,t0,x0,t_end,h_first,tolerance);

% Print and save the initial point and all accepted points.
filename = 'ch09_rkf45_scalar.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open data file: %s',filename);
end
fprintf('Adaptive Runge-Kutta-Fehlberg RKF45\n');
fprintf('Tolerance: %.4e\n',tolerance);
fprintf('%10s %18s %14s\n','t','numerical x','h');
fprintf(fid,'# t x h\n');
for i = 1:length(t)
    fprintf('%10.6f %18.12f %14.6e\n',t(i),x(i),h(i));
    fprintf(fid,'%.16e %.16e %.16e\n',t(i),x(i),h(i));
end
fclose(fid);
fprintf('Accepted steps: %8d\n',accepted);
fprintf('Rejected trials:%8d\n',rejected);
fprintf('Final x:       %18.12f\n',x(end));
if check_exact
    exact = atan(100*t_end); % DEFAULT example only
    fprintf('Exact final x: %18.12f\n',exact);
    fprintf('Abs error:     %12.4e\n',abs(x(end)-exact));
end
fprintf('Data file: %s\n',filename);

% Solution: dense exact curve and actual numerical points.
figure;
if check_exact
    t_plot = linspace(t0,t_end,2000);
    x_plot = atan(100*t_plot);
    plot(t_plot,x_plot,'k-','LineWidth',1.2);
    hold on;
    plot(t,x,'ko--','MarkerSize',4);
    hold off;
    legend('Exact','RKF45','Location','best');
else
    plot(t,x,'ko-','MarkerSize',4);
end
xlabel('t');
ylabel('x(t)');
title('Adaptive RKF45: solution');
grid on;

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

%==================== Equation and numerical method =====================
function dx = f(t,x)
% Supply x' = f(t,x) here. x is available for state-dependent equations.
dx = 100/(1+10000*t*t);
end

function [t_values,x_values,h_values,accepted,rejected] = ...
    rkf45(f,t0,x0,t_end,h_first,tolerance)
% Adaptive Fehlberg RKF45. Advance with the fifth-order estimate.
% The fourth/fifth-order difference estimates local error; tolerance
% is not a bound on the accumulated global error.
% A rejected trial leaves both time and solution unchanged.
% Returned column vectors contain the initial and accepted points only.
max_attempts = 100000;
h_min = 1.0e-12;
if t_end <= t0 || h_first <= 0 || tolerance <= 0
    error('Use t_end > t0, h_first > 0, and tolerance > 0.');
end
t = t0;
x = x0;
h = h_first;
accepted = 0;
rejected = 0;
attempts = 0;
t_values = t0;
x_values = x0;
h_values = 0;

while t < t_end
    % Ignore only a machine-rounding gap at the endpoint.
    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;

    % Fehlberg stages: each k includes h.
    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;
    err = abs(k1/360-128*k3/4275-2197*k4/75240+k5/50+2*k6/55);
    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_values(end+1,1) = x;
        h_values(end+1,1) = h;
    else
        rejected = rejected+1;
    end

    % Embedded error scales approximately as h^5.
    % Safety factor 0.9; limit growth and shrinkage.
    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