I would suggest down-sampling. Make a run an find estimates, then focus on those orbits of interest with more data points from your set.

I am using artificial directed networks and first I test on some small size network, run the Kuramoto model on that and calculate the LEs to get some insights.

I did not get how to focus on some orbits if interest, could you please introduce some reference?

You mean just using time series or using the Jocobian matrix?

I feel the purpose of using large n network should be addressed. But, you can follow Abolfazl idea on small size network and compared the trend as you increase the size. Check the uploaded papers if they are useful.

Thank you Adedayo Oke Adelakun

Large network is required to study emergent behaviors.

The first one was helpful.

A possible alternative is to consider evolution of the "overall" distance between neighboring solutions instead of multidimensional perturbation vector. It can be done in the following way:

1) take some initial condition "of interest";

2) impose a small random perturbation of initial conditions on each node, and calculate the initial "overall" distance as D(0) = sqrt(sum_j |d_j|^2);

3) integrate equations of motion for unperturbed and perturbed multidimensional trajectories over small (!) time interval Delta t, and find D(Delta t);

4) calculate time-local Lyapunov exponent as lambda_i = (1/Delta t)log(D(Delta t)/D(0));

5) assign the final state of the unperturbed trajectory as new unperturbed initial condition;

6) repeat steps 2-5 many-many times (number of repetitions depends on your particular interests) and then find the "global" Lyapunov exponent by averaging lambda_i.

A similar procedure was used in Article Chaos-assisted formation of immiscible matter-wave solitons ...

Thank you @Denis for your answer.

It would be nice to have a rough estimate of the time complexity of the method.

Time complexity is nearly the same as for integration of two individual trajectories. All other calculations are not very time-consuming.