Solving the system purely numerically (without any symbolic computations) is a bit faster. This can be done by using the signed incidence matrix of the graph (matrix B below).

I have to use the \[Theta]_?VectorQ pattern in otder to prevent NDSolve from starting any symbolic analysis.

ParametricNDSolveValue provides a neat way to collect everyting into a single function.

Also, there is no point in solving the ODE until tmax if the solution is evaluated only at evaltime; so I let just return the state at that time, not the whole InterpolationFunction.

ClearAll[\[Lambda], T];
B = N[IncidenceMatrix[Graph[UpperTriangularize[A]["NonzeroPositions"],DirectedEdges -> True]]];
X[\[Theta]_?VectorQ, \[Lambda]_] := N[omegadata] - \[Lambda] B . Sin[\[Theta] . B];
newa = ParametricNDSolveValue[{
  \[Theta]'[t] == X[\[Theta][t], \[Lambda]], 
  \[Theta][0] == \[Theta]start}, \[Theta][T], 
  {t, 0, T}, {\[Lambda], T}];

aa = a[0.2, evaltime]; // AbsoluteTiming // First
bb = newa[0.2, evaltime]; // AbsoluteTiming // First
Max[Abs[aa - bb]/Max[Abs[aa]]]

1.65085

0.019422

4.63889*10^-8

This turns out to be almost 100 times faster. Probably, a different ODE integrator is applied, which explains the small difference of the results.

Solving the system purely numerically (without any symbolic computations) is a bit faster. This can be done by using the signed incidence matrix of the graph (matrix B below).

I have to use the \[Theta]_?VectorQ pattern in otder to prevent NDSolve from starting any symbolic analysis.

ParametricNDSolveValue provides a neat way to collect everyting into a single function.

Also, there is no point in solving the ODE until tmax if the solution is evaluated only at evaltime; so I let just return the state at that time, not the whole InterpolationFunction.

ClearAll[\[Lambda], T];
B = N[IncidenceMatrix[Graph[UpperTriangularize[A]["NonzeroPositions"],DirectedEdges -> True]]];
X[\[Theta]_?VectorQ, \[Lambda]_] := N[omegadata] - \[Lambda] B . Sin[\[Theta] . B];
newa = ParametricNDSolveValue[{
  \[Theta]'[t] == X[\[Theta][t], \[Lambda]], 
  \[Theta][0] == \[Theta]start}, \[Theta][T], 
  {t, 0, T}, {\[Lambda], T}];

aa = a[0.2, evaltime]; // AbsoluteTiming // First
bb = newa[0.2, evaltime]; // AbsoluteTiming // First
Max[Abs[aa - bb]/Max[Abs[aa]]]

1.65085

0.019422

4.63889*10^-8

This turns out to be almost 100 times faster. Probably a different ODE integrator is applied, which explains the small difference of the results.