This is a really good question. Any attempt at a good answer requires breaking the news that in a realistic solar system with some big planets "A Keplerian orbit is merely an idealized, mathematical approximation at a particular time". If you are writing a computer program and not using ink, feathers, and parchment, it's best to just integrate the appropriate ODEs for the six element state vectors ($x, y, z, v_x, v_y, v_z$) of each body in question.
There are several ways to do this. You can use an ephemeris to get the state vectors at one moment in time (set up initial conditions) and integrate yourself, or you can get tables of positions of everything besides your asteroid and interpolate them. Two examples of resources for these are NASA JPL's Horizons and the python package Skyfield but there are several other options.
You can do some looking around here and in the Space Exploration stackexchange sites for tags for those.
The sun's motion relative to the solar system's center of mass is not a simple orbit, but reflects the motion of (primarily) Jupiter, Saturn, Uranus and Neptune, and the only way to get it right is to either calculate all of them yourself, or use a pre-calculated ephemeris.
These are necessary both to get the position of the sun correctly, and to calculate the gravitational effects of those planets on your asteroid!
Remember though, that gravity travels at the speed of light, so if you want to do precise calculations you have to at least make a first order approximation at each time step to where the planet would have been in the past such that the gravity "signal" will have just arrived at your asteroid at the time reflected in your time step.
Here's an example of the sun's motion with python. Dashed lines represent the equatorial radius of the sun - the motion regularly moves beyond one solar radius from the solar system center of mass.
import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader
load = Loader('~/Documents/fishing/SkyData')
data = load('de405.bsp') # this one is big, de421 is smaller but more limited years
ts = load.timescale()
years = np.arange(1980, 2030, 0.1, dtype=float)
times = ts.utc(years, 1, 1)
sun = data['sun']
sunpos = sun.at(times).ecliptic_position().km
sunlims = 695700. * np.array([-1, 1])[:, None] * np.ones_like(years)
if 1 == 1:
plt.figure()
for thing in sunpos:
plt.plot(years, 1E-06*thing)
plt.ylabel('millions of km', fontsize=16)
plt.xlabel('year', fontsize=16)
for sunlim in sunlims:
plt.plot(years, 1E-06*sunlim, '--k')
plt.title('J2000 ecliptic position of the Sun (x, y, z)', fontsize=16)
plt.show()
"""
# you can save as a numpy or json file:
np.save('sunpos', sunpos) # You can save as .npy if you want to
with open('sunpos.json', 'w') as outfile:
json.dump(sunpos.round(1).tolist(), outfile)
"""