If you're not doing N-body physics, so the planets do not interact then you can use analytic solutions to the Kepler problem. Eventually you'll realize that you need to solve hyperbolic orbits at some point as well. That will lead you to universal variables formulations of solving the Kepler problem.

The best solutions to that are probably going to be Goodyear's method:

W. Goodyear, “Completely General Closed Form Solution for Coordinates and Partial Derivatives of the Two-Body Problem”, The Astronomical Journal, Vol. 70, No. 3, 1965, pp. 189–192 (or the NASA NTRS TD document on the same material)

Shepperd's method:

Shepperd, S.W. Celestial Mechanics (1985) 35: 129. https://doi.org/10.1007/BF01227666

Or Danby-Stumpff:

Danby, J.M.A. Celestial Mechanics (1987) 40: 303. https://doi.org/10.1007/BF01235847

There is some MATLAB code here which might be useful (and vastly more accessible), although random code snippets on matlabcentral are far from guaranteed to be bug free and it looks like this code may lack useful normalization of its inputs (generally you're going to want to normalize to the scale of your problem so that you do math in units where r0-bar = 1.0 and mu-bar = 1.0 and where v-bar = 1 is the velocity in a circular orbit at r0 or something like that).

If you are going to do N-body integration of planetary motion then I think you're going to have to use numerical integration. Runge-Kutta will violate conservation of Energy so you will likely want to use Symplectic Integration. The 4th order symplectic integrator in that article is not that difficult to code -- although that leaves you with the difficulty of guessing the correct timestep (again, normalization helps because a circular planetary orbit and circular LEO are the same problem just with different distance scales) and with interpolation of the interior points (and you need to watch out for Runge's phenomenon, but I haven't wrestled with that, so don't know which approach to take there).

If you're going to use Runge-Kutta then Dormand-Prince with dynamic step side and its 3rd order interpolant will be very convenient, and is what Matlab uses in its ode45 solver.

I would probably advise starting with the simplest runge-kutta implementation based on ease of coding, but if you're doing runge-kutta on every physics tick to advance it forward one step then that is pretty brutal and the errors will eventually add up, but you could prototype it that way. At some point you'll want to go to a system where you solve the problem for many time steps into the future, and then you use an interpolating function to pick off the solution at intermediate timesteps (which is the point of my mentioning Dormand-Prince and its interpolating function).

While there has already been a high quality accepted answer for years, here is some additional background, some particularly helpful resources, and additional tips for first-time orbit propagation.

If you're not doing N-body physics, so the planets do not interact then you can use analytic solutions to the Kepler problem. Eventually you'll realize that you need to solve hyperbolic orbits at some point as well. That will lead you to universal variables formulations of solving the Kepler problem.

The best solutions to that are probably going to be Goodyear's method:

W. Goodyear, “Completely General Closed Form Solution for Coordinates and Partial Derivatives of the Two-Body Problem”, The Astronomical Journal, Vol. 70, No. 3, 1965, pp. 189–192 (or the NASA NTRS TD document on the same material)

Shepperd's method:

Shepperd, S.W. Celestial Mechanics (1985) 35: 129. https://doi.org/10.1007/BF01227666

Or Danby-Stumpff:

Danby, J.M.A. Celestial Mechanics (1987) 40: 303. https://doi.org/10.1007/BF01235847

There is some MATLAB code here which might be useful (and vastly more accessible), although random code snippets on matlabcentral are far from guaranteed to be bug free and it looks like this code may lack useful normalization of its inputs (generally you're going to want to normalize to the scale of your problem so that you do math in units where r0-bar = 1.0 and mu-bar = 1.0 and where v-bar = 1 is the velocity in a circular orbit at r0 or something like that).

If you are going to do N-body integration of planetary motion then I think you're going to have to use numerical integration. Runge-Kutta will violate conservation of Energy so you will likely want to use Symplectic Integration. The 4th order symplectic integrator in that article is not that difficult to code -- although that leaves you with the difficulty of guessing the correct timestep (again, normalization helps because a circular planetary orbit and circular LEO are the same problem just with different distance scales) and with interpolation of the interior points (and you need to watch out for Runge's phenomenon, but I haven't wrestled with that, so don't know which approach to take there).

If you're going to use Runge-Kutta then Dormand-Prince with dynamic step side and its 3rd order interpolant will be very convenient, and is what Matlab uses in its ode45 solver.

I would probably advise starting with the simplest runge-kutta implementation based on ease of coding, but if you're doing runge-kutta on every physics tick to advance it forward one step then that is pretty brutal and the errors will eventually add up, but you could prototype it that way. At some point you'll want to go to a system where you solve the problem for many time steps into the future, and then you use an interpolating function to pick off the solution at intermediate timesteps (which is the point of my mentioning Dormand-Prince and its interpolating function).