In this paper, we study the following nonlinear boundary value problem D \(\begin{aligned} \left\{ \begin{array}{ll} D\psi -a(x)\psi =f(x,\psi )+\epsilon h(x,\psi )&{} \quad \hbox {on }M \\ B_{CHI}\psi =0 &{} \quad \hbox {on }\partial M \end{array} \right. \end{aligned}\) where M is a compact Riemannian spin manifold of dimension \(m\ge 2\) and the boundary \(\partial M\) has non-negative mean curvature, and D is the Dirac-Atiyah-Singer operator. Let S(M) denote the spinor bundle on M and \(\psi :M\rightarrow S(M)\) be a section. a(x) is a scalar field on M, \(\epsilon \in \mathbb {R}\) . \(f(x,\psi )\) , \(h(x,\psi ):M\times S_{m}\rightarrow S_{m}\) are two nonlinear map, and \(B_{CHI}\) is the chirality boundary operator. Under some mild assumptions on a, f, and h, we obtain the ground state solution of (D) with \(\epsilon =0\) and infinitely many large norm solutions of (D) with \(\epsilon \in \mathbb {R}\) and infinitely many small energy solutions of (D) with \(\epsilon >0\) by using variational methods.