close all; clear; clc;

%%
b = 0.1; h = 1/2; n = 3;
X = exp(1i*(0:(n-1))*2*pi/n);

S = @(x) sum((X-x(1))./(abs(X-x(1)).^2+h^2).^(5/2));
dxdt = @(t,x) [0 1;-1 -b]*x + [0;1]*S(x(1));

rs = linspace(-2,2,2^5); is = rs;
C = zeros(numel(is),numel(rs)); % Preallocate
as = randperm(numel(C)); % Random starting points
for ai = 1:numel(as)
    [ri,ii] = ind2sub(size(C),as(ai)); % Convert to indices
    x0 = [rs(ri)+is(ii).*1i; 0];
    [t,y] = ode23(dxdt,[0 20],x0);
    [~,C(ii,ri)] = min(abs(X-y(end,1)));
    if ~mod(ai,10)
        figure(1); clf;
        subplot(1,2,1); cla; hold on;
        plot(y(:,1)); plot(y(1,1),'o');  plot(y(end,1),'p'); plot(X,'k.');
        axis square; box on; grid on;
        xlim([-1 1]*2); ylim([-1 1]*2);
        subplot(1,2,2); cla;
        imagesc(C);
        colormap(inferno);
        axis tight square; grid on;
        drawnow;
        
        fprintf('%.2f%% done\n',ai/numel(as)*100);
    end
end
