We propose an iterative method to safely learn the unmodeled dynamics of a nonlinear system using Bayesian Gaussian process (GP) models with polynomial kernel functions. The method maintains safety by ensuring that the system state stays within the region of attraction (ROA) of a stabilizing control policy while collecting data. A quadratic programming based exploration control policy is computed to keep the exploration trajectory inside an inner-approximation of the ROA and to maximize the information gained from the trajectory. A prior GP model, which incorporates prior information about the unknown dynamics, is used to construct an initial stabilizing policy. As the GP model is updated with data, it is used to synthesize a new policy and a larger ROA, which increases the range of safe exploration. The use of polynomial kernels allows us to compute ROA inner-approximations and stabilizing control laws for the model using sum-of-squares programming. We also provide a probabilistic guarantee of safety which ensures that the policy computed using the learned model stabilizes the true dynamics with high confidence.