  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}=3\lambda\phi^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.
