  In this work, we present an analytical solution for QCD$\otimes$QED coupled
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations at the
leading order (LO) accuracy in QED and next-to-leading order (NLO) accuracy in
perturbative QCD using double Laplace transform. This technique is applied to
obtain the singlet, gluon and photon distribution functions and also the proton
structure function. We also obtain contribution of photon in proton at LO and
NLO at high energy and successfully compare the proton structure function with
HERA data \cite{12} and APFEL results \cite{7}. Some comparisons also have been
done for the singlet and gluon distribution functions with the MSTW
results\cite{9}. In addition, the contribution of photon distribution function
inside the proton has been compared with results of MRST \cite{11} and with the
contribution of sea quark distribution functions which obtained by MSTW
\cite{9} and CTEQ6M \cite{14}.
