Why do objects in a gravity simulation experience sudden large accelerations?
If you really want a general gravitation simulator (i.e. one that will handle more then two bodies), then there are methods for reducing the error involved in the simulation, but there aren't any methods for eliminating the error. Below are a few approaches - none of these approaches are perfect, since there's a balance between physical accuracy, programming & computational speed, and reducing the escaping behavior.
Adaptive time steps - when your particle is experiencing large changes in acceleration (i.e. when it's moving fast near a massive body), you can reduce the time step $dt$ so you're "skipping" less time. This will increase the physical accuracy of the simulation, but it will also slow down the simulation (and slow it down unevenly), so it's not a good approach if you want it to look good in real time.
Hard spheres - make your masses hard spheres that bounce off each other. This will reduce the buildup of error, since the masses won't get as close to each other. This can increase the physical accuracy of the simulation, unless you object to planets that bounce like rubber balls (which is a fair objection, I suppose - you could also make them inelastic spheres, which is probably more accurate for planets).
Speed limit - you could just program in a hard speed limit to limit escaping behavior. Whenever the velocity exceeds the speed limit, resize the velocity magnitude back to the speed limit. This isn't very physically accurate, and could result in some strange looking behavior, but it's easy to program, and will reduce the escaping of your masses.
Conservation of energy - at each time step, calculate the total amount of gravitational and kinetic energy all the masses have. At each time step, if the total amount of energy has changed, artificially adjust the velocities of the masses so the total amount of energy stays the same. This is also not perfectly accurate, but it does maintain fidelity with one physical law, and it will reduce the escaping behavior.
If you'd like help understanding the implementation of one of these methods, I can explain in more detail.
Brionus has touched on the key - adaptive time steps. When you start getting large accelerations, reduce the size of your time increments. Also, when you are not accelerating much, increase the size.
A fairly standard way to do this is to calculate your position change over one step. Then cut the step in half and, starting from the same starting point, compute 2 successive position steps. Compare the two final position changes, and if they differ by some predetermined factor (let's say 10^-6), replace the original time step with the smaller step, and do the calculation all over again. If the two steps matched closely, try a computation with a time step twice the original.
This takes a lot of extra computation, but it produces a simulation which is neither too precise nor not precise enough. For orbital simulations, the large time steps used during low-gravity trajectories will more than compensate for the extra computation time.
EDIT - In response to the request for "the calculus approach":
The calculation approach which you used will work for a really simple simulator, but it ignores the interactions between the non-sun bodies, as well as assuming that the central sun does not move. To handle this, use a more generalized approach. Store the (x,y,z) (Vx,Vy,Vz) and (Mx,My,Mz) values for your bodies in arrays indexed the same way. Then the gravitational attractions between any two bodies will simply be calculated as (Fx,Fy,Fz), where $$F_n =\frac{\Delta n}{\sqrt{{\Delta x}^2 + {\Delta y}^2 + {\Delta z}^2}} \frac{GM_a M_b}{{\Delta x}^2 + {\Delta y}^2 + {\Delta z}^2}$$ Calculate each component separately, and integrate separately to get new velocities and positions. Also note that you will need to calculate $N(N-1)$ values for N bodies (including the sun) but that body A pulls on body B exactly as hard as body B pulls on body A, so you only need to do the numerical computation half of the apparent total, although you'll need to be careful to reverse the sign to get the other half of the values right. (The reason it's $N(N-1)$ rather than $N^2$ is that you don't calculate the pull of a body on itself.)
For a simple simulator, you can use Euler integration. Given the force on a body and its mass, for a very small $\Delta t$ you can say $$\Delta V = \frac{F \Delta t}{M}$$ and when you know both the original velocity V and the change in velocity, the change in position $\Delta P$ is$$\Delta P = {(V + \Delta V)\Delta t}$$
For shortish time scales and very small $\Delta t$ this will work, but you are fundamentally trying to approximate an irregular curve with straight segments and in the long run you'll see large and growing errors. So for longer periods you'll want to get into more sophisticated algorithms. As zeldridge has commented, Runge-Kutta is a well-known alternative.
SECOND EDIT - Also, when updating your values, do each set of calculations based on the same original conditions for all bodies. That is, if you're calculating Vnew and Pnew for body A, do not calculate the results for body B using the updated Pnew for body A. Calculate an entire new array of V and P values, then replace the old values as a block.
If you are not interested in a full n-body simulation and will accept one where the sun is much more massive than all the planets, you can simplify things a lot. We have an analytic solution to the two body problem, so you can apply the sun's gravity that way. For each planet, given its position and velocity, you can compute its angular momentum and energy, then get the orbit. This lets you compute the location of the planet at the end of the time step ignoring the attractions of all the planets. To zeroth order, you are done, because you ignore the planet-planet interactions. If you want to apply the planet-planet interactions, you can compute the force on planet $i$ from all the other planets and apply the acceleration for the entire time step. As the force is much smaller than the sun's gravity, the change in position will be small. This will keep your errors much smaller, roughly as the ratio of the forces on the planet from the other planets to the force from the sun. You can then apply the adaptive time steps and energy conservation that the other responders suggest.