Prove that $f(x, y) \le 3 $ for $x \ge 0, y > 0$ 
Let $x \ge 0, y>0$ and
\begin{align*}
f(x,y)&=\sqrt{\dfrac{y}{y+x^2}}+4\sqrt{\dfrac{y}{(y+(x+1)^2)(y+(x+3)^2)}}\\[6pt]
&\qquad +4\sqrt{\dfrac{y}{(y+(x-1)^2)(y+(x-3)^2)}}.
\end{align*}
Prove that $f(x,y) \le 3$.

I can prove when $x\ge 2, f(x,y) < 3$, but $f(x,y)=3$ when $x=0,y=3$, so the key part is $x\le 2 $
I found when $x<2, f'_x(y \ge 1.5) <0$ but I can't prove it as it has  high degree equations.
And for $x<2,y<1.5$, I have no idea how to prove $f(x,y)<3$

No clue without an overview. With help of restrictions as found already by the OP, a contour plot / isoline chart has been produced for the function at hand:

The $y$-axis is in $\color{green}{\mbox{green}}$. Our viewport is:

  xmin := -6 ; xmax := 6;
  ymin := 0 ; ymax := 12;

The are 27 contour levels nivo, defined by:


  for geval := 1 to 27 do
  begin
    nivo := geval/9.1;

The higher the level, the more black. The lower the level, the more white. 
The pixels in $\color{red}{red}$ are close to the conjectured maximum $=3$ . They are defined by:


const
  eps : double = 0.0002;
  if (3-f < eps) and (3-f > 0) then

It is seen that the red spot is inside an area with very black lines,
meaning that the function $f(x,y)$ is increasing there, towards the red spot.

Analysis. There is a mirror symmetry in the $y$-axis, as is clear from the figure but also from some algebra.
Thus for all real $x$ and for all real $y > 0$ we have:
$$
f(x,y) = f(-x,y)
$$
It follows that:
$$
\frac{\partial f(x,y)}{\partial x} = - \frac{\partial f(-x,y)}{\partial x}
\quad \Longrightarrow \quad
\left.\frac{\partial f(x,y)}{\partial x}\right|_{x=0} = - \left.\frac{\partial f(-x,y)}{\partial x}\right|_{x=0}
\\ \Longrightarrow \quad
\left.\frac{\partial f(x,y)}{\partial x}\right|_{x=0} = 0
$$
Meaning that all function gradients near the $y$-axis are tangent to it. Since the gradient vectors are perpendicular
to the isolines, this can be observed in the figure as well. Now a sufficient condition for having a stationary point
(maximum or minimum eventually) is that the total gradient be zero. But we already have $\,\partial f / \partial x = 0\,$
at the $y$-axis. Therefore substitute $\,x=0\,$ in $\,f(x,y)\,$ and only consider function values at the $y$-axis:
$$
f(0,y)=1+8\sqrt{\frac{y}{(y+1)(y+9)}}
$$
Extremes are found for:
$$
\frac{d f(0,y)}{dy} = 0 \quad \Longrightarrow \\
\frac{1}{2} \left[ {\frac {1}{ \left( y+1 \right)  \left( y+9 \right) }}-{
\frac {y}{ \left( y+1 \right) ^{2} \left( y+9 \right) }}-{\frac {y}{
 \left( y+1 \right)  \left( y+9 \right) ^{2}}} \right] 
\sqrt{\frac{\left( y+1 \right)  \left( y+9 \right)}{y}} = 0 \\
\Longrightarrow \quad (y+1)(y+9) - y(y+9) - y(y+1) = -y^2 +9 = 0 \quad \Longrightarrow \quad y = \pm 3
$$
With the restriction $\,y > 0$ . Hence the only extreme at $\,(y,f(0,y))$ is $(3,3)\,$ and it is a maximum.
Picking up the comment by nbubis, prove that $f(0,2) < 3$ and $f(0,4) < 3$:
$$
f(0,2)=1+8\sqrt{\frac{2}{(2+1)(2+9)}}=1+\sqrt{\frac{128}{33}}<1+\sqrt{\frac{128}{32}}\\
f(0,4)=1+8\sqrt{\frac{4}{(4+1)(4+9)}}=1+\sqrt{\frac{256}{65}}<1+\sqrt{\frac{256}{64}}
$$
