We establish a uniform error estimate of a finite difference method for the Klein-Gordon-Schrodinger (KGS) equations with two dimensionless parameters \begin{document}$0 and \begin{document}$0 , which are the mass ratio and inversely proportional to the speed of light, respectively. In the simultaneously nonrelativistic and massless limit regimes, i.e., \begin{document}$γ\sim\varepsilon$\end{document} and \begin{document}$\varepsilon \to 0^+$\end{document} , the KGS equations converge singularly to the Schrodinger-Yukawa (SY) equations. When \begin{document}$0 , due to the perturbation of the wave operator and/or the incompatibility of the initial data, which is described by two parameters \begin{document}$α≥0$\end{document} and \begin{document}$β≥-1$\end{document} , the solution of the KGS equations oscillates in time with \begin{document}$O(\varepsilon)$\end{document} -wavelength, which requires harsh meshing strategy for classical numerical methods. We propose a uniformly accurate method based on two key points: (ⅰ) reformulating KGS system into an asymptotic consistent formulation, and (ⅱ) applying an integral approximation of the oscillatory term. Using the energy method and the limiting equation via the SY equations with an oscillatory potential, we establish two independent error bounds at \begin{document}$O(h^2+τ^2/\varepsilon)$\end{document} and \begin{document}$O(h^2+τ^2+τ\varepsilon^{α^*}+\varepsilon^{1+α^*})$\end{document} with \begin{document}$h$\end{document} mesh size, \begin{document}$τ$\end{document} time step and \begin{document}$α^* = \min\{1, α, 1+β\}$\end{document} . This implies that the method converges uniformly and optimally with quadratic convergence rate in space and uniformly in time at \begin{document}$O(τ^{4/3})$\end{document} and \begin{document}$O(τ^{1+\frac{α^*}{2+α^*}})$\end{document} for well-prepared ( \begin{document}$α^* = 1$\end{document} ) and ill-prepared ( \begin{document}$0≤α^* ) initial data, respectively. Thus the \begin{document}$\varepsilon$\end{document} -scalability of the method is \begin{document}$τ = O(1)$\end{document} and \begin{document}$h = O(1)$\end{document} for \begin{document}$0 , which is significantly better than classical methods. Numerical results are reported to confirm our error bounds. Finally, the method is applied to study the convergence rates of KGS equations to its limiting models in the simultaneously nonrelativistic and massless limit regimes.