%% Animated trajectories with moving particles for the ODE
%   x' = v
%   v' = (1 - 3x^2 - 2v^2)*v - x

clear; close all; clc;

%% Domain
xmin = -1.5; xmax = 1.5;
ymin = -1.5; ymax = 1.5;

%% Figure setup
figure('units','normalized','outerposition',[0 0 1 1]);
hold on;
axis equal;
xlabel('x'); ylabel('v');
title('Animated Trajectories of the ODE System');
xlim([xmin xmax]); ylim([ymin ymax]);
set(gca,'FontSize',12,'Box','on');

%% ODE system
odefun = @(t,z) [ ...
    z(2);                                   % x' = v
    (1 - 3*z(1)^2 - 2*z(2)^2)*z(2) - z(1)    % v'
];

%% Initial conditions on boundary of square
N = 5;
b = linspace(xmin, xmax, N)';

IC_boundary = [
    b , ymin*ones(N,1);     % bottom
    b , ymax*ones(N,1);     % top
    xmin*ones(N,1) , b;     % left
    xmax*ones(N,1) , b      % right
];

%% === Initial conditions near the origin ===
% Option A: small circle around origin
theta = linspace(0, 2*pi, 8)';
r = 0.01;
IC_circle = [ r*cos(theta), r*sin(theta) ];

% Option B: small grid near origin
[xg, vg] = meshgrid(linspace(-0.2,0.2,3), linspace(-0.2,0.2,3));
IC_grid = [ xg(:), vg(:) ];

% Choose one (or combine):
IC_near_origin = [IC_circle; IC_grid];

%% Combine all ICs
IC = [IC_boundary; IC_near_origin];

num_traj = size(IC,1);

%% Time settings
tspan = [0 30];
t_anim = linspace(tspan(1), tspan(2), 1000);

Xsol = zeros(num_traj, length(t_anim));
Vsol = zeros(num_traj, length(t_anim));

%% Solve ODE trajectories
for k = 1:num_traj
    z0 = IC(k,:);
    [t,z] = ode45(odefun, tspan, z0);
    Xsol(k,:) = interp1(t, z(:,1), t_anim, 'linear', 'extrap');
    Vsol(k,:) = interp1(t, z(:,2), t_anim, 'linear', 'extrap');
end

%% Create graphics objects
trajLines = gobjects(num_traj,1);
dots = gobjects(num_traj,1);

colors = lines(num_traj);

for k = 1:num_traj
    trajLines(k) = plot(nan, nan, 'Color', colors(k,:), 'LineWidth', 1.4);
    dots(k) = plot(Xsol(k,1), Vsol(k,1), 'o', ...
        'MarkerFaceColor', colors(k,:), 'Color', colors(k,:), 'MarkerSize', 10);
end

%% Animation loop
disp('Press any key to start animation...');
pause;

for j = 1:length(t_anim)
    if mod(j,3) == 1
        for k = 1:num_traj
            set(trajLines(k), 'XData', Xsol(k,1:j), ...
                              'YData', Vsol(k,1:j));
            set(dots(k), 'XData', Xsol(k,j), ...
                         'YData', Vsol(k,j));
        end
        drawnow;
    end
end