%=======================================================================
% Chapter 10, code 0: shooting for a second-order boundary value problem.
% Translation of ch10_shooting.f90, with solution plots.
% Solve y'' = f(x,y,y'), with y(a) = alpha and y(b) = beta.
% Secant iteration finds the unknown initial slope. Each shooting trial
% is an initial value problem integrated with adaptive Fehlberg RKF45.
%
% Example: y'' = y, y(0) = y(1) = 1.
% Exact solution: cosh(x-0.5)/cosh(0.5).
% MATLAB alternatives: ode45 can integrate each shooting trial (it uses
% Dormand-Prince, rather than Fehlberg RKF45). bvp4c solves the BVP directly
% by collocation and requires an initial solution guess.
%
% To solve another problem:
%   1. Edit the interval, boundary values, and slope guesses below.
%   2. Edit f(x,y,dy), immediately below the driver.
%   3. Disable check_exact or replace exact_solution.
% Returned Y has one grid point per row: Y(:,1) = y, Y(:,2) = y'.
%=======================================================================
clear;

%-------------------- Problem and controls: edit here --------------------
a = 0.0;
b = 1.0;
alpha = 1.0;                 % y(a)
beta = 1.0;                  % y(b)
slope1 = 0.0;                % First guess for y'(a)
slope2 = 1.0;                % Second guess for y'(a)
h_first = 0.05;
h_max = 0.05;                % Maximum spacing between output points
ivp_tolerance = 1.0e-10;      % Absolute LOCAL tolerance for RKF45
boundary_tolerance = 1.0e-8;  % Required |y(b)-beta|
capacity = 100001;           % Maximum stored points per shooting trial
check_exact = true;
make_plot = true;
%------------------------------------------------------------------------
fprintf('Boundary value problem: shooting with secant iteration and RKF45\n');
fprintf('RKF45 tolerance:   %12.4e\n',ivp_tolerance);
fprintf('Boundary tolerance:%12.4e\n',boundary_tolerance);
[x,Y,slope,residual,trials,accepted,rejected] = shoot2(@f,a,b,alpha,beta, ...
    slope1,slope2,h_first,h_max,ivp_tolerance,boundary_tolerance,capacity);

% Print and save only the converged trajectory, including both endpoints.
filename = 'ch10_shooting.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open data file: %s',filename);
end
fprintf('\n');
if check_exact
    fprintf('%10s %16s %16s %16s %12s\n', ...
        'x','numerical y','dy/dx','exact y','abs error');
    fprintf(fid,'# x y dy_dx exact_y abs_error\n');
else
    fprintf('%10s %16s %16s\n','x','numerical y','dy/dx');
    fprintf(fid,'# x y dy_dx\n');
end
max_error = 0.0;
for i = 1:length(x)
    if check_exact
        exact = exact_solution(x(i));
        abs_error = abs(Y(i,1)-exact);
        max_error = max(max_error,abs_error);
        fprintf('%10.6f %16.8e %16.8e %16.8e %12.4e\n', ...
            x(i),Y(i,1),Y(i,2),exact,abs_error);
        fprintf(fid,'%24.16e %24.16e %24.16e %24.16e %24.16e\n', ...
            x(i),Y(i,1),Y(i,2),exact,abs_error);
    else
        fprintf('%10.6f %16.8e %16.8e\n',x(i),Y(i,1),Y(i,2));
        fprintf(fid,'%24.16e %24.16e %24.16e\n',x(i),Y(i,1),Y(i,2));
    end
end
fclose(fid);
fprintf('\nShooting trials:      %8d\n',trials);
fprintf('Initial slope:        %16.10f\n',slope);
fprintf('Boundary residual:    %16.6e\n',residual);
fprintf('Accepted steps:       %8d\n',accepted);
fprintf('Rejected trials:      %8d\n',rejected);
if check_exact
    fprintf('Largest sampled error:%16.6e\n',max_error);
end
fprintf('Step counts refer to the final shooting trial.\n');
fprintf('Data file: %s\n',filename);

% Plot the solution and its derivative. A dense exact curve is used only
% for comparison; numerical markers are the actual accepted grid points.
if make_plot
    figure;
    subplot(2,1,1);
    plot(x,Y(:,1),'ko-','MarkerSize',3,'LineWidth',1.0);
    if check_exact
        hold on;
        x_plot = linspace(a,b,500)';
        plot(x_plot,exact_solution(x_plot),'k--','LineWidth',1.2);
        hold off;
        legend('Numerical','Exact','Location','best');
    end
    xlabel('x');
    ylabel('y(x)');
    title('Boundary value problem: shooting solution');
    grid on;
    subplot(2,1,2);
    plot(x,Y(:,2),'ko-','MarkerSize',3,'LineWidth',1.0);
    xlabel('x');
    ylabel('dy/dx');
    title('Derivative of the solution');
    grid on;
end


function value = f(x,y,dy)
% The original second-order equation: y'' = f(x,y,y').
% x and dy are available for equations that depend on them.
value = y;
end


function value = exact_solution(x)
% Analytic solution for the supplied example and boundary values.
% Accepts a scalar or a vector, including the dense plotting grid.
value = cosh(x-0.5)/cosh(0.5);
end


function [grid,Y,slope,residual,trials,accepted,rejected] = shoot2(fcn,a,b, ...
    alpha,beta,slope1,slope2,h_first,h_max,ivp_tolerance,boundary_tolerance,capacity)
%=======================================================================
% Treat s = y'(a) as an unknown. Each trial starts with y(a) = alpha,
% y'(a) = s. The mismatch R(s) = y(b;s)-beta must vanish.
% Two successive slopes and residuals determine the next secant slope.
% The final integration returns the converged trajectory.
% Secant iteration is not bracketed; convergence depends on the guesses.
% A small boundary residual alone does not establish solution accuracy:
% tighten the RKF45 tolerance and compare trajectories for a new problem.
%=======================================================================
max_trials = 30;
if boundary_tolerance <= 0 || slope1 == slope2
    error('Use positive boundary tolerance and two distinct slope guesses.');
end
fprintf('%8s %18s %18s\n','trial','initial slope','boundary residual');
slope = slope1;
previous_slope = slope1;
previous_residual = 0.0;
for trials = 1:max_trials
    initial = [alpha; slope];
    [grid,Y,final,accepted,rejected] = rkf45_system(fcn,a,initial,b, ...
        h_first,h_max,ivp_tolerance,capacity);
    residual = final(1)-beta;
    fprintf('%8d %18.10e %18.10e\n',trials,slope,residual);
    if abs(residual) <= boundary_tolerance
        return;
    end
    if trials == max_trials
        break;
    end
    if trials == 1
        next_slope = slope2;
    else
        denominator = residual-previous_residual;
        scale = max([abs(residual),abs(previous_residual),realmin]);
        if abs(denominator) <= eps*scale
            error('Secant residuals are indistinguishable; change the guesses.');
        end
        next_slope = slope-residual*(slope-previous_slope)/denominator;
        if ~isfinite(next_slope)
            error('Nonfinite secant slope.');
        end
        if next_slope == slope
            error('Secant slope cannot advance; check tolerances.');
        end
    end
    previous_slope = slope;
    previous_residual = residual;
    slope = next_slope;
end
error('Shooting did not converge within the trial limit.');
end


function dy = equations(fcn,x,y)
% Two first-order equations: y(1)' = y(2), y(2)' = f(x,y(1),y(2)).
dy = [y(2); fcn(x,y(1),y(2))];
if ~all(isfinite(dy))
    error('Nonfinite derivative in trial integration.');
end
end


function [grid,Y,xf,accepted,rejected] = rkf45_system(fcn,t0,x0,t_end, ...
    h_first,h_max,tolerance,capacity)
%=======================================================================
% Adaptive Fehlberg RKF45. All step decisions are inside this routine.
% Save initial and accepted points; trial trajectories produce no output.
% Advance with the fifth-order estimate. The largest absolute embedded
% component error controls acceptance. One absolute tolerance assumes
% components of comparable magnitudes and controls LOCAL error.
% Internal states are two-element columns; output Y has one point per row.
%=======================================================================
max_attempts = 100000;
if capacity < 2 || capacity ~= floor(capacity) || t_end <= t0 ...
        || h_first <= 0 || h_max <= 0 || tolerance <= 0
    error('Use integer capacity >= 2, b > a, positive steps, and positive tolerance.');
end
grid = zeros(capacity,1);
Y = zeros(capacity,2);
t = t0;
h = h_first;
state = x0;
accepted = 0;
rejected = 0;
attempts = 0;
grid(1) = t;
Y(1,:) = state.';
while t < t_end
    if attempts >= max_attempts
        error('Maximum number of trials exceeded.');
    end
    h = min([h,h_max,t_end-t]);
    tf = t+h;
    if tf <= t
        error('Step too small to advance the independent variable.');
    end
    attempts = attempts+1;
    % Every k includes h. Both equations use the SAME intermediate state.
    k1 = h*equations(fcn,t,state);
    work = state+k1/4;
    k2 = h*equations(fcn,t+h/4,work);
    work = state+3*k1/32+9*k2/32;
    k3 = h*equations(fcn,t+3*h/8,work);
    work = state+1932*k1/2197-7200*k2/2197+7296*k3/2197;
    k4 = h*equations(fcn,t+12*h/13,work);
    work = state+439*k1/216-8*k2+3680*k3/513-845*k4/4104;
    k5 = h*equations(fcn,t+h,work);
    work = state-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40;
    k6 = h*equations(fcn,t+h/2,work);
    trial = state+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);
    if ~all(isfinite(trial)) || ~all(isfinite(component_error))
        error('Nonfinite RKF45 result: check the model and slope guesses.');
    end
    err = max(component_error);
    if err <= tolerance
        if h == t_end-t
            tf = t_end;
        end
        t = tf;
        state = trial;
        accepted = accepted+1;
        if accepted+1 > capacity
            error('Increase trajectory storage capacity.');
        end
        grid(accepted+1) = t;
        Y(accepted+1,:) = state.';
    else
        rejected = rejected+1;   % Keep both t and state unchanged.
    end
    if err == 0
        factor = 2.0;
    else
        factor = 0.9*(tolerance/err)^0.2;
    end
    if err <= tolerance
        h = h*min(2.0,factor);
    else
        h = h*max(0.2,factor);
    end
end
xf = state;
grid = grid(1:accepted+1);
Y = Y(1:accepted+1,:);
end