A new technique for computing the reachable set of hybrid systems with nonlinear continuous dynamics is presented. Previous work showed that abstracting the nonlinear continuous dynamics to linear differential inclusions results in a scalable approach for reachability analysis. However, when the abstraction becomes inaccurate, linearization techniques require splitting of reachable sets, resulting in an exponential growth of required linearizations. In this work, the nonlinearity of the dynamics is more accurately abstracted to polynomial difference inclusions. As a consequence, it is no longer guaranteed that reachable sets of consecutive time steps are mapped to convex sets as typically used in previous works. Thus, a non-convex set representation is developed in order to better capture the nonlinear dynamics, requiring no or much less splitting. The new approach has polynomial complexity with respect to the number of continuous state variables when splitting can be avoided and is thus promising when a linearization technique requires splitting for the same problem. The benefits are presented by numerical examples.