function plot_tracer_migration()
% =========================================================================
% TRACER MIGRATION IN SQUEEZE FLOW (analytical; no fitted parameters)
% Integrates the transverse tracer dynamics
%     dzeta/ds  = zeta(1-zeta)(1-2zeta),   s = ln(h0/h),  a = exp(s/2)
%     dln(r)/ds = 3 zeta (1-zeta)
% and reproduces the tracer-kinematics figure: (a) the local rate profile
% H(zeta)/Hbar = 6 zeta (1-zeta), (b) convergence of zeta onto the
% mid-plane attractor, (c) radial trajectories versus the mid-plane law
% r ~ a^(3/2) and the depth-averaged law r ~ a.
%
% Every ode45 trajectory is cross-checked against the exact closed form
% (Sect. "Exact solution and the attractor structure"):
%     eps(s) = eps0 / sqrt(4 eps0^2 + (1-4 eps0^2) e^s)
%     r(s)/r0 = [4 eps0^2 + (1-4 eps0^2) e^s]^(3/4)
% and the asymptotic memory factor M(zeta0) = [4 zeta0 (1-zeta0)]^(3/4)
% is reported (r/r0 -> M(zeta0) a^(3/2) as s -> infinity).
% =========================================================================

S     = 4.0;                                  % total compression, a = e^{S/2}
z0s   = [0.05 0.15 0.30 0.50 0.70 0.85 0.95]; % initial relative heights
sgrid = linspace(0, S, 400);
opts  = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);

figure('Color', 'w', 'Position', [60, 120, 1350, 420]);

% ---- (a) local expansion-rate profile ----------------------------------
subplot(1, 3, 1); hold on;
zz = linspace(0, 1, 300);
plot(6*zz.*(1-zz), zz, 'k-', 'LineWidth', 2);
plot([1 1], [0 1], '--', 'Color', [0 0.45 0.74], 'LineWidth', 1.5);
plot([0 1.75], [0.5 0.5], ':', 'Color', [0.84 0.15 0.16], 'LineWidth', 1.5);
plot(1.5, 0.5, 'o', 'MarkerFaceColor', [0.84 0.15 0.16], ...
     'MarkerEdgeColor', 'k', 'MarkerSize', 8);
text(1.02, 0.90, 'gap average $\bar H$', 'Interpreter', 'latex', ...
     'Color', [0 0.45 0.74], 'FontSize', 9, 'Rotation', 90);
text(1.55, 0.53, 'mid-plane: $1.5\bar H$', 'Interpreter', 'latex', ...
     'Color', [0.84 0.15 0.16], 'FontSize', 9);
xlabel('$H(\zeta)/\bar{H} = 6\zeta(1-\zeta)$', 'Interpreter', 'latex');
ylabel('$\zeta = z/h$', 'Interpreter', 'latex');
title('(a) rate across the gap', 'Interpreter', 'latex');
xlim([0 1.75]); ylim([0 1]); grid on;

% ---- (b),(c) trajectory integration, cross-checked against closed form --
cols = parula(numel(z0s));
hb = gobjects(1, numel(z0s));   % line handles, panel (b), for the legend
hc = gobjects(1, numel(z0s));   % line handles, panel (c), for the legend
leg_labels = arrayfun(@(z) sprintf('$\\zeta_0=%.2f$', z), z0s, ...
                       'UniformOutput', false);

subplot(1, 3, 2); hold on;
subplot(1, 3, 3); hold on;
fprintf('  zeta0    zeta(S)     r/r0    max|exact-ode45|    M(zeta0)\n');
maxerr_all = 0;
for i = 1:numel(z0s)
    [~, y] = ode45(@rhs, sgrid, [z0s(i); 0], opts);
    av = exp(sgrid/2);

    % exact closed form, Eq. (exactsol)
    e0  = z0s(i) - 0.5;
    D   = 4*e0^2 + (1 - 4*e0^2) .* exp(sgrid);
    zex = 0.5 + e0 ./ sqrt(D);
    rex = D .^ 0.75;
    err = max(abs([zex(:) - y(:,1), rex(:) - exp(y(:,2))]), [], 'all');
    maxerr_all = max(maxerr_all, err);

    subplot(1, 3, 2);
    hb(i) = plot(av, y(:, 1), '-', 'Color', cols(i, :), 'LineWidth', 1.6);
    subplot(1, 3, 3);
    hc(i) = plot(av, exp(y(:, 2)), '-', 'Color', cols(i, :), 'LineWidth', 1.6);

    Mz = (4*z0s(i)*(1 - z0s(i)))^0.75;   % asymptotic memory factor
    fprintf('  %5.2f    %6.3f   %8.3f    %10.2e         %7.4f\n', ...
            z0s(i), y(end, 1), exp(y(end, 2)), err, Mz);
end
fprintf('  depth-averaged a       = %8.3f\n', exp(S/2));
fprintf('  mid-plane law a^(3/2)  = %8.3f\n', exp(S/2)^1.5);
fprintf('  max |exact - ode45| over all trajectories and s: %.3e\n', maxerr_all);

% ---- (b) mid-plane attractor --------------------------------------------
subplot(1, 3, 2);
plot([1 exp(S/2)], [0.5 0.5], ':', 'Color', [0.84 0.15 0.16], 'LineWidth', 1.8);
set(gca, 'XScale', 'log'); xlim([1 exp(S/2)]); grid on;
xlabel('$a = \sqrt{h_0/h}$', 'Interpreter', 'latex');
ylabel('$\zeta(a)$', 'Interpreter', 'latex');
title('(b) mid-plane attractor', 'Interpreter', 'latex');
legend(hb, leg_labels, 'Interpreter', 'latex', 'Location', 'eastoutside', ...
       'FontSize', 8);

% ---- (c) radial trajectories ---------------------------------------------
subplot(1, 3, 3);
av = exp(sgrid/2);
hmid   = plot(av, av.^1.5, '--', 'Color', [0.84 0.15 0.16], 'LineWidth', 2.2);
hfield = plot(av, av,      '-.', 'Color', [0 0.45 0.74],   'LineWidth', 2.0);
set(gca, 'XScale', 'log', 'YScale', 'log'); xlim([1 exp(S/2)]); grid on;
xlabel('$a = \sqrt{h_0/h}$', 'Interpreter', 'latex');
ylabel('$r(a)/r_0$', 'Interpreter', 'latex');
title('(c) radial trajectories', 'Interpreter', 'latex');
legend([hc, hmid, hfield], [leg_labels, ...
       {'$r \propto a^{3/2}$ (mid-plane)', '$r \propto a$ (field)'}], ...
       'Interpreter', 'latex', 'Location', 'eastoutside', 'FontSize', 8);

    function dy = rhs(~, y)
        zt = y(1);
        dy = [zt*(1 - zt)*(1 - 2*zt); 3*zt*(1 - zt)];
    end
end