N_pendulum_multi_data

clc;clear;close all
load('func_N3k.mat')
opts = odeset('RelTol',1e-13,'AbsTol',1e-16,'MaxStep',5e-3, ...
    'InitialStep',1e-5,'Stats','on');
T=100;fps=30;Frame=T*fps;
t=linspace(0,T,Frame);
M=100;
x=zeros(Frame,N+1,M);y=x;
tic
for ii=1:M
u0=[pi/2+zeros(N,1)-ii/M*1e-8;zeros(N,1)];
[~,u] = ode89(@(t,u) eqns(t,u,N,f),t,u0,opts);
for k=1:N
    x(:,k+1,ii)=x(:,k,ii)+sin(u(:,k));
    y(:,k+1,ii)=y(:,k,ii)-cos(u(:,k));
end
end
toc
function du=eqns(~,u,N,f)
du=[u(N+1:2*N);
    f(u(1:N))*([1;u(2*N:-1:N+1).^2])];
end