For the system of Fredholm integral equations \[ u_i(t)= \int_0^1g_i(t,s)f_i\left(s,u_1(s),u_2(s),\dots,u_n(s)\right)\,ds,\quad t\in[0,1],\quad 1\leq i\leq n,\tag{1} \] the authors provide conditions on the nonlinearities \(f_i\) and kernels \(g_i\), so that the system (1) has at least three positive solutions. The main tool is the Leggett-Williams fixed-point theorem. The obtained results are applied to the system of second order Neumann boundary value problems \[ \begin{aligned} -u''_i(t)+c_iu_i(t)&=f_i(t,u_1(t),u_2(t),\dots,u_n(t)),\quad t\in[0,1],\\ u_i'(0)&=u_i'(1)=0,\quad 1\leq i\leq n,\\ \end{aligned} \] where \(c_i>0\) are fixed. Similarly, the results can be also applied to other well-known systems of boundary value problems, for example, \((m,p)\), Lidstone, focal, conjugate, Hermite and Sturm-Liouville, and extend the corresponding work in the literature.