  We analytically solved the QED $\otimes$ QCD coupled DGLAP evolution
equations at leading order (LO) quantum electrodynamics (QED) and next to
leading order (NLO) quantum chromodynamics (QCD) approximations, using the
Laplace transform method and then computed the proton structure function in
terms of the unpolarized parton distributions functions. Our analyitical
solutions for parton densities are in good agreement with those from APFEL (A
PDF Evolution Library) (Computer Physics Communications 185, 1647-1668 (2014))
and CT14QED (Phys. Rev. D 93, 114015 (2016)) global parameterizations. We also
compared the proton structure function, $F_{2}^{p}(x,Q^{2})$, with experimental
data released by the ZEUS and H1 collaborations at HERA. There is a nice
agreement between them in the range of low and high x and $Q^{2}$.
