lorenz_system.hpp

  using state_type = std::array<long double, 3>;
  // Lorenz attractor
  auto lorenz_system = [&,p = 10.L,r = 28.L,b = 8.L/3](const long double t, const state_type& x) {
    state_type dxdt{};
    // x=x[0], y=x[1], z=x[2]
    dxdt[0] = -p*x[0] + p*x[1];              // dx/dt
    dxdt[1] = -x[0] * x[2] + r*x[0] - x[1];  // dy/dt
    dxdt[2] = x[0] * x[1] - b*x[2];          // dz/dt
    return dxdt;
  }