This paper examines a bias correction method to decrease inaccuracies of SGP4 propagation for cubesats. Using only historical TLE data and TLE differencing methods, residuals of SGP4 propagation are fit into first and second order polynomials and exponential equations for various coordinate components from TEME, RSW, NTW, PQW coordinate systems as well as classical orbital, equinoctial and flight elements. A total of 36 different elements are used. The most appropriate type of fit equation for each test case is selected by comparing root mean squared errors. A bias detection concept based on confidence intervals is developed to programmatically understand if error behaviour is systematic (biased) in a certain coordinate element. It is shown that certain coordinate elements suffer systematic errors most of the time while error growth on others are seemingly random, although complete behaviour is different for each satellite. Estimations of error growth equations for coordinate elements that suffer systematic errors are used to modify SGP4 results to reduce propagation errors. Success, bias detection and error reduction rates of this bias correction operation are investigated with respect to propagation duration, fit span and goodness of fit. It has been found that drastically increasing the accuracy of SGP4 propagation (for both position and velocity predictions) is possible with efficient selection of coordinate elements and statistical properties. As a result, it is shown that for certain coordinate elements, if a systematic error growth is detected, results of SGP4 can be modified and modifications lead to improvements for more than 99% of the biased cases. It is shown that improvements can be small (less than 1%) or very large (up to 99%) with a mean of 52% error reduction for all coordinate elements.