We analyze reachability properties and local input/output gains of systems with polynomial vector fields. Upper bounds for the reachable set and nonlinear system gains are characterized by polynomial Lyapunov (storage) functions satisfying certain bilinear constraints. A methodology utilizing information from simulations to generate Lyapunov function candidates satisfying necessary conditions for bilinear constraints is proposed. The suitability of Lyapunov function candidates are assessed solving linear sum-of-squares optimization problems. Qualified candidates are used to compute upper bounds for the reachable set and nonlinear system gains and to initialize further coordinate-wise affine optimization. We illustrate the method on several examples from the literature.