  We construct spherically-symmetric static solutions of the
Einstein-Klein-Gordon-Euler system involving a complex scalar field governed by
a periodic potential which emerges in models of axion-like particles, and
fermionic matter modeled by a perfect fluid with a polytropic equation of
state. Such solutions describe gravitationally bound composites of fermions and
axions which we dub as fermion-axion stars. Sequences of pure axion-stars in
the existence domain may show the presence of multiple stable branches
depending on the value of the decay constant parameter in the potential; this
reflects in the appearance of multiple islands of stability in the
2-dimensional parameter space of fermion-axion configurations. We investigate
the domain of existence for three different values of the decay constant,
identifying one or more regions of linear stability making use of a method we
already employed in previous works. We confirm the results from the linear
analysis performing fully non-linear numerical relativity evolutions. In this
context we perform several numerical simulations to identify regions in the
parameter space where unstable models face different fates: the collapse to a
Schwarzschild black hole, the migration to a stable model and finally the
dispersion of the scalar field together with the dilution of the fermionic
matter, which approaches a static fermion star model with very low mass. This
latter scenario was never observed in previous models without the periodic
potential.
