In $f(R)$ gravity and Brans-Dicke theory with scalar potentials, we study the structure of neutron stars on a spherically symmetric and static background for two equations of state: SLy and FPS. In massless BD theory, the presence of a scalar coupling $Q$ with matter works to change the star radius in comparison to General Relativity, while the maximum allowed mass of neutron stars is hardly modified for both SLy and FPS equations of state. In Brans-Dicke theory with the massive potential $V(phi)=m^2 phi^2/2$, where $m^2$ is a positive constant, we show the difficulty of realizing neutron star solutions with a stable field profile due to the existence of an exponentially growing mode outside the star. As in $f(R)$ gravity with the $R^2$ term, this property is related to the requirement of extra boundary conditions of the field at the surface of star. For the self-coupling potential $V(phi)=lambda phi^4/4$, this problem can be circumvented by the fact that the second derivative $V_{,phi phi}=3lambdaphi^2$ approaches 0 at spatial infinity. In this case, we numerically show the existence of neutron star solutions for both SLy and FPS equations of state and discuss how the mass-radius relation is modified as compared to General Relativity.