% fast_slow_animation_corrected.m
% System:
%   xdot = (1/eps)*( y - F(x) )
%   ydot = -eps * x
% With F(x) = x^3 - x

clear; close all; clc;

%% Parameters
eps = .05;                 % epsilon
F = @(x) x.^3 - x;          % nonlinearity F(x)

tmax = 200;                 % long simulation time
dt_out = 0.2;              % output time step
slowdown = 8;               % animation slowdown factor

%% Initial condition
x0 = 1.5;
y0 = 1.5;
ic = [x0; y0];

%% Integrator
opts = odeset('RelTol',1e-6, 'AbsTol',1e-9);
odefun = @(t,z) [ (1/eps)*( z(2) - F(z(1)) );  -eps * z(1) ];

tspan = 0:dt_out:tmax;
[~,Z] = ode15s(odefun, tspan, ic, opts);
x = Z(:,1);  
y = Z(:,2);
t = tspan(:);

%% Figure setup
figure('units', 'normalized', 'outerposition', [0,0,1,1]);

%% ----------------- Phase Portrait -----------------
ax1 = subplot(1,2,1);
hold(ax1,'on'); box(ax1,'on');
xlabel(ax1,'x'); ylabel(ax1,'y');
title(ax1,'Phase plane');

% x-nullcline: y = F(x)
xgrid = linspace(-3, 3, 400);
plot(ax1, xgrid, F(xgrid), 'k-', 'LineWidth',1.6);

% y-nullcline: x = 0
plot(ax1, [0 0], [-3 3], 'r--', 'LineWidth',1.6);

legend(ax1, 'x-nullcline: y=F(x)', 'y-nullcline: x=0', 'Location','northwest');

% Phase space limits
xlim(ax1, [-3 3]); 
ylim(ax1, [-3 3]);

% Full trajectory reference
plot(ax1, x, y, 'Color',[0.8 0.8 0.8]);

% animated trail + moving dot
trail_len = round(1/dt_out);
hTrail = plot(ax1, NaN, NaN, 'b-', 'LineWidth',1.8);
hDot   = plot(ax1, x(1), y(1), 'ro', 'MarkerFaceColor','r', 'MarkerSize',8);

%% ----------------- x(t) Panel -----------------
ax2 = subplot(1,2,2);
hold(ax2,'on'); box(ax2,'on');
xlabel(ax2,'t'); ylabel(ax2,'x(t)');
title(ax2,'x(t)');

xlim(ax2, [t(1), t(end)]);
ylim(ax2, [min(x)-0.5, max(x)+0.5]);

% faint full trace
plot(ax2, t, x, 'Color',[0.85 0.85 0.85]);

% animated x(t) trace
hXt   = plot(ax2, t(1), x(1), 'b-', 'LineWidth',1.8);
hXdot = plot(ax2, t(1), x(1), 'ro', 'MarkerFaceColor','r', 'MarkerSize',8);

drawnow;

%% ----------------- Animation Loop -----------------
nSteps = length(t);
tic;
for k = 1:nSteps
    
    % phase portrait trail
    idx0 = max(1, k-trail_len);
    set(hTrail, 'XData', x(idx0:k), 'YData', y(idx0:k));
    set(hDot,   'XData', x(k),      'YData', y(k));
    
    % x(t)
    set(hXt,   'XData', t(1:k), 'YData', x(1:k));
    set(hXdot, 'XData', t(k),   'YData', x(k));

    % titles
    title(ax1, sprintf('Phase plane (t = %.2f)', t(k)));
    title(ax2, sprintf('x(t) (t = %.2f)', t(k)));

    drawnow limitrate;

    % slowdown  
    elapsed = toc;
    target = t(k)/slowdown;
    if target > elapsed
        pause(min(0.05, target - elapsed));
    end
end

fprintf('Animation finished at t = %.2f.\n', tmax);