The strong force has an approximate symmetry called isospin, which is a rotation in an abstract space spanned by the neutron and proton. This symmetry is broken by the Coulomb force, which only acts on protons, and by the slight mass difference between up quarks and down quarks, which makes the neutron slightly heavier than the proton.
Consider possible bound states of neutrons and protons. These are fermions, so the wave function has to be anti-symmetric. Neutrons and protons are heavy, so we can analyze the situation using non-relativistic quantum mechanics. If a bound state exists, the ground state is in an s-wave (symmetric) because the centrifugal barrier is purely repulsive. This means that the spin-isospin wave function has to be anti-symmetric. This leaves two possibilities:
$$
(I=0,S=1) \;\;or\;\; (I=1,S=0)
$$
where $I$ and $S$ are the total isopsin and spin of the state. Note that $I,S=0$ is anti-symmetric in isospin/spin, and $I,S=1$ is symmetric. The $I=0$ state corresponds to
$$
|I=0\rangle = \frac{1}{\sqrt{2}}\left(|np\rangle -|pn\rangle\right)
$$
and $I=1$ has three isospin components
$$
|I=1,I_3=+1,0,-1\rangle = \left\{|pp\rangle ,\frac{1}{\sqrt{2}}\left(|np\rangle +|pn\rangle\right),|nn\rangle\right\}
$$
The potential in the $I=0,1$ channels is determined by complicated interactions between quarks and gluons.
Empirically (or from numerical calculations in lattice QCD) we know that both the $I=0$ potential and the $I=1$ potential are attractive, but the $I=0$ interaction is slightly more attractive.
Note that in non-relativistic QM (in three dimensions), an attractive interaction is not sufficient to form a bound state. The potential has to have a minimum strength. Nuclear physics is fine-tuned, the empirical NN potentials are close to the threshold where bound states appear.
There is an $I=0$ bound state, called the deuteron (binding energy 2.2 MeV), but the di-neutron (and its isospin partners) fall just short of being bound. Due to the Coulomb force, there is a little extra repulsion in $pp$, but the $nn$ channel is almost bound.
This can be quantified using the scattering length. The $nn$ scattering length is about 20 fm, and the corresponding energy scale is
$$
B=\frac{1}{2ma^2}\simeq 50 \, keV
$$
so the di-neutron comes within 50 keV of being bound. This also means that if one could change the quark masses slightly, there would presumably be a bound di-neutron (this can be checked in lattice QCD).
It is possible in principle to have a three neutron bound state even if there is no two neutron bound state (states like that are called Borromean), but this is disfavored by the Pauli principle (not all three neutrons can be in an s-wave). The same applies to 4n bound states (there is a recent claim of a 4n resonance http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.116.052501).
In practice $^3H=pnn$ and $^3He=ppn$ exist, but not $3n$.
Note that neutron stars, clusters of $10^{57}$ neutrons, do exist. Of course, these objects are gravitationally bound.
I suspect you are relying on the modern language, which is yet controverted by the effective theory community these days, if I am not too cut off from recent developments... I believe it is all hiding behind the receding obsession with renormalizability, and thus minimal coupling, obviated by the Wilson revolution.
The point is the minimal-coupling gauge-invariant, renormalizable, Dirac action was perfectly adequate to describe g=2 through the Gordon decomposition term of the current associated with the magnetic dipole density of the electron,
$$
-j'_\mu A^\mu\sim -(e/2m) \left (\frac{1}{2} F^{\mu\nu} \bar \psi \sigma _{\mu\nu}\psi\right ).
$$
A physicist of the late 30s (I'm cluelessly guessing here!), knowing the nucleon magnetic moments were not canonical, would augment his minimal coupling Dirac action for them with an extra, nonminimal-coupling (unrenormalizable, which he wouldn't know about) Pauli moment term, stuck in by hand,
$$
-(e/2M) \left (\frac{1}{2} F^{\mu\nu} \bar \psi \sigma _{\mu\nu}\psi\right ),
$$
perhaps to be added to the above Gordon current piece (which would vanish for the neutral neutron! whose magnetic moment was measured by Alvarez & Bloch, 1939), for a phenomenological parameter M. He would fit everything to determine M for experimentally determined nuclear magnetons; note is not the nucleon's mass, but merely of its rough order of magnitude; and hope for the future to clarify things. Not having a clue about the mysteries of nature, he'd leave it at that.
The late 40s revolution in renormalization allowed computation of corrections to the g of the electron; but, due to non-renormalizability, not for the nucleon, what with the above nasty dimension-5 spatchcocked Pauli term, with its mystery scale M. (As an aside, this term is dear to the heart of extended-supergravitists, M being the Planck scale.)
Then, in the mid-60s, during the triumphant march of quark composition, the said Pauli moment terms were further calculated from a loosely-bound constituent-quark wavefunction. I would not be surprised if today's lattice geeks can specify the exact parameters in the effective Dirac action-cum-Pauli moment involved.
In the subsequent SM years, launched by 't Hooft's proof of SSB YM renormalizability, giddy quasi-religious attachment to renormalizability was lavished on these systems—until Ken Wilson restored humility by reminding us we all live in a resolutely effective action world. But, "elementary" was a virtual shorthand for a field described by a renormalizable action.
So, in the early 80s, hyperambitious model builders were ready to contemplate compositeness even for pure Dirac action particles like the leptons, included in the sample papers of my comment above, and Harari 1982. Now they had the reverse problem: how to constrain the scales of compositeness, so in effect, how to make the M of an extraneous Pauli moment be enormous. I hope you are not asking about that, since these guys went pretty deeply pretty fast. And then they appeared to mumble, shrug, and walk away.
Best Answer
The one-word answer is yes.
You are also correct that the neutron is not just a proton and electron living together. The process of merging a proton and electron proceeds via the weak force. Specifically, an up quark in the proton exchanges a W boson with the electron. The W boson carries a unit of positive charge from the quark to the electron. In that process the up quark (charge +2/3) is converted to a down quark (charge -1/3) so that the proton (uud) becomes a neutron (udd). The negatively charged electron is converted into a neutrino. This is one important point left out in your question. The full reaction is $p+e^-\to n+\nu_e$.
There is a general principle of quantum field theory called crossing symmetry that roughly states that for any process I can exchange what I call initial and final particles. So you are correct that neutron decay $n\to p+ e^- + \bar\nu_e$ implies that the process $p+e^-\to n+\nu_e$ can also happen.
This process does also happen in nature. It is one mode of radioactive decay of nuclei. Some nuclei with a sufficiently large number of protons can become more stable by absorbing one of their electrons and converting one proton into a neutron. This can happen because electron orbitals have a small but non-zero overlap with the nucleus, so that they "sometimes come into contact with" the protons.
This process can also happen artificially as you suggest. In fact, it seems that accelerators used in medical facilities produce neutrons as a by-product, exactly as you suggest, and this is apparently a difficulty that must be dealt with, see this paper.
In general, because the mass difference between the proton and neutron is about an MeV, in any system including protons and electrons at a temperature of order an MeV or higher, there will necessarily be populations of both neutrons and protons connected to each other by such processes, with relative amounts determined by the relevant Boltzmann factors. This should include systems where thermal fusion is taking place.
However the actual process of producing helium from hydrogen, as far as I understand, does not depend on capturing an electron on a proton to form a neutron. In stellar nucleosynthesis, two protons merge to form deuterium. That is, in the process of merging, one proton is converted into a neutron by the emission of a positron and a neutrino. Helium-2 (two protons) is highly unstable, so that this proton-to-neutron conversion producing stable deuterium is more important.