SGP4 is the standard procedure that TLEs are intended to work with. All of the following are extremely helpful, but the most important point would be use a standard, recent SGP4 package that is recommended, do not try to use the elements in a TLE yourself. This is becuase the TLE and the SGP4 package are built for each other.
One interesting point is that for high altitude orbits with periods longer than 225 minutes, a modern SGP4 implementation will switch to an algorithm that is called SDP4. See the question “Deep space” corrections in SGP4; how does it account for the Sun's and Moon's gravity? for more on that.
The easiest to access SGP4 that I know is in the Python package Skyfield. You can find SGP4 routines in many languages in many places. I would recommend you choose something that is used widely and well maintained, and not just any random code on the internet.
Skyfield's SGP4 is based on :https://github.com/brandon-rhodes/python-sgp4
Here's a plot of the low earth orbit of the SpaceX F9 upper stage known as Roadster generated using Skyfield and Roadster's TLE, of course before the second burn that put it in heliocentric orbit.

def makecubelimits(axis, centers=None, hw=None):
lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
if centers == None:
centers = [0.5*sum(pair) for pair in lims]
if hw == None:
widths = [pair[1] - pair[0] for pair in lims]
hw = 0.5*max(widths)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)
print("hw was None so set to:", hw)
else:
try:
hwx, hwy, hwz = hw
print("ok hw requested: ", hwx, hwy, hwz)
ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
except:
print("nope hw requested: ", hw)
ax.set_xlim(centers[0]-hw, centers[0]+hw)
ax.set_ylim(centers[1]-hw, centers[1]+hw)
ax.set_zlim(centers[2]-hw, centers[2]+hw)
return centers, hw
TLE = """1 43205U 18017A 18038.05572532 +.00020608 -51169-6 +11058-3 0 9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()
from skyfield.api import Loader, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180
load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts = load.timescale()
planets = load('de421.bsp')
earth = planets['earth']
Roadster = EarthSatellite(L1, L2)
print(Roadster.epoch.tt)
hours = np.arange(0, 3, 0.01)
time = ts.utc(2018, 2, 7, hours)
Rpos = Roadster.at(time).position.km
Rposecl = Roadster.at(time).ecliptic_position().km
print(Rpos.shape)
re = 6378.
theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
cph, sph = [f(phi) for f in (np.cos, np.sin)]
lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
lon0[1]*cph + lon0[0]*sph,
lon0[2]) )
lons.append(lon)
lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
cph, sph = [f(phi) for f in (np.cos, np.sin)]
lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
lats.append(lat)
if True:
fig = plt.figure(figsize=[10, 8]) # [12, 10]
ax = fig.add_subplot(1, 1, 1, projection='3d')
x, y, z = Rpos
ax.plot(x, y, z)
for x, y, z in lons:
ax.plot(x, y, z, '-k')
for x, y, z in lats:
ax.plot(x, y, z, '-k')
centers, hw = makecubelimits(ax)
print("centers are: ", centers)
print("hw is: ", hw)
plt.show()
r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re
if True:
plt.figure()
plt.plot(hours, r_Roadster)
plt.plot(hours, alt_roadster)
plt.xlabel('hours', fontsize=14)
plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
plt.show()