Improve performance of is_energy_preserving via sparse matrices

[ranocha]

This tackles the biggest chunk of #130.

Current main:

julia> using BSeries

julia> begin # second run after compilation

       # rational coefficients
       series = bseries(AverageVectorFieldMethod(), 10)
       @time is_energy_preserving(series)

       # floating point coefficients
       series = bseries(AverageVectorFieldMethod(Float64), 10)
       @time is_energy_preserving(series)

       end
  0.664049 seconds (119.29 k allocations: 88.295 MiB, 0.62% gc time)
  0.560968 seconds (118.06 k allocations: 80.759 MiB, 0.35% gc time)
true

This PR:

julia> using BSeries

julia> begin # second run after compilation

       # rational coefficients
       series = bseries(AverageVectorFieldMethod(), 10)
       @time is_energy_preserving(series)

       # floating point coefficients
       series = bseries(AverageVectorFieldMethod(Float64), 10)
       @time is_energy_preserving(series)

       end
  0.397417 seconds (111.10 k allocations: 13.948 MiB)
  0.327424 seconds (109.87 k allocations: 14.197 MiB)
true

With this PR, modified_equation takes 93% of the total runtime of is_energy_preserving(series) for series = bseries(AverageVectorFieldMethod(Float64), 10).