In this paper, we study the existence of positive solutions to the boundary value problem for the fractional differential system $$\left\{\begin{array}{lll} D_{0^+}^\beta \phi_p(D_{0^+}^\alpha u) (t) = f_1 (t, u (t), v (t)),\quad t \in (0, 1),\\ D_{0^+}^\beta \phi_p(D_{0^+}^\alpha v) (t) = f_2 (t, u (t), v(t)), \quad t \in (0, 1),\\ D_{0^+}^\alpha u(0)= D_{0^+}^\alpha u(1)=0,\; u (0) = 0, \quad u (1)-\Sigma_{i=1}^{m-2} a_{1i}\;u(\xi_{1i})=\lambda_1,\\ D_{0^+}^\alpha v(0)= D_{0^+}^\alpha v(1)=0,\; v (0) = 0, \quad v (1)-\Sigma_{i=1}^{m-2} a_{2i}\; v(\xi_{2i})=\lambda_2, \end{array}\right. $$ where \({1<\alpha,\beta\leq 2, 2 <\alpha + \beta\leq 4, D_{0^+}^\alpha}\) is the Riemann–Liouville fractional derivative of order α. By using the Leggett–Williams fixed point theorem in a cone, the existence of three positive solutions for nonlinear singular boundary value problems is obtained.