close all; clear; clc

pendulums = 200;

X = zeros(2 * pendulums, 1);
X(1:2:end) = linspace(-2.9, 2.9, pendulums);

A = diag(ones(1, 2 * pendulums - 1), 1);
A(2:2:end, :) = 0;

B = diag(ones(1, 2 * pendulums - 1), -1);
B = B(:, 1:2:end);

dt = 0.001;
t = 0:dt:10;

g = 9.8;
l = 1;

plot_data = cell(1, pendulums);

figure

for j = 1:pendulums
    plot_data{j} = animatedline();
end

xlim([0, 10]);
ylim([-3, 3]);

for i = 1:length(t)
    X = X + dt * (A * X - (g / l) * B * sin(X(1:2:end)));
    if mod(i, 50) == 0
        for j = 1:pendulums
            addpoints(plot_data{j}, t(i), X(2 * j - 1));
        end
        drawnow limitrate
    end
end