An analytical solution based on the Laplace transformation technique for the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi DGLAP evolution equations is presented at next-to-leading order accuracy in perturbative QCD. This technique is also applied to extract the analytical solution for the proton structure function, F p 2 (x, Q 2 ), in the Laplace s-space. We present the results for the separate parton distributions for all parton species, including valence quark densities, the anti-quark and strange sea parton distribution functions (PDFs), and the gluon distribution. We successfully compare the obtained parton distribution functions and the proton structure function with the results from GJR08 [Eur. Phys. J C 53 (2008) 355-366] and KKT12 [J. Phys. G 40 (2013) 045002] parametrization models as well as the x-space results using QCDnum code. Our calculations show a very good agreement with the available theoretical models as well as the deep inelastic scattering (DIS) experimental data throughout the small and large values of x. The use of our analytical solution to extract the parton densities and the proton structure function is discussed in detail to justify the analysis method considering the accuracy and speed of calculations. Overall, the accuracy we obtain from the analytical solution using the inverse Laplace transform technique is found to be better than 1 part in 10 4 to 10 5 . We also present a detailed QCD analysis of non-singlet structure functions using all available DIS data to perform global QCD fits. In this regard we employ the Jacobi polynomial approach to convert the results from Laplace s space to Bjorken x space. The extracted valence quark densities are also presented and compared to the JR14, MMHT14, NNPDF and CJ15 PDFs sets. We evaluate the numerical effects of target mass corrections (TMCs) and higher twist (HT) terms on various structure functions, and compare fits to data with and without these corrections.