Conservation of energy is universal across all mechanics, including orbital mechanics. The orbital energy will be constant as long as there are no non-conservative external forces acting on the body (gravity is a conservative force). Semi-major axis is directly related to orbital energy as you have shown, therefore it will also remain constant. The semi-major axis defines the size of the orbit and therefore it also makes sense that the orbit doesn't change without any external forces.
If you are only working on a two-body problem, then you can simply choose an orbit based on the set of orbital elements (semi-major axis, eccentricity, and argument of perigee, as well as inclination and right ascension of the ascending node in full 3D cases). With those elements you don't have to simulate the orbit with gravity forces and equations of motion, you can just simply map out the position of the body at any time you wish.
However, let's assume you still want to simulate gravity for the two-body case since it can be useful and is easily extended to more than two bodies. Looking at your code, the biggest thing that stands out to me is: why are you dividing the distance by 3000? Remove that scaling from your function and only apply it for plotting.