First of all, let us show that the function $f$ is monotone increasing in the wide sense, namely it satisfies the property that $a < b \implies f(a) \le f(b)$. So, suppose $a < b$. Then we can choose a number $t > 0$ such that $b = t^2 + a + t f(4a)$. Substituting $x = t$, $y = a$ into the given functional equation, we obtain $f(t^2) + f(a) = f(b)$. Since $f(t^2) \ge 0$, we then obtain $f(a) \le f(b)$, showing the desired monotone property of $f$. Furthermore, we note that the same argument shows that if $f(t) > 0$ for an arbitrary $t > 0$, then $a < b \implies f(a) < f(b)$.

We separate the following two cases:
(1): $f(t) > 0$ is satisfied for an arbitrary $t > 0$.
(2): $f(t) = 0$ for some $t > 0$.

Let us first determine $f(x)$ in the case (1). As we stated above, in this case, we have $a < b \implies f(a) < f(b)$, and therefore, we must have $f(a) = f(b) \implies a = b$. ($f$ is said to be injective if it satisfies this property.) For arbitrary non-negative real numbers $a$ and $b$, substitute $(x, y) = (\sqrt{a}, b)$, $(\sqrt{b}, a)$ into the given functional equation, then we obtain
$$\begin{align*} f(a) + f(b) &= f(a + b + \sqrt{a} f(4b)), \\ f(b) + f(a) &= f(b + a + \sqrt{b} f(4a)) \end{align*}$$
from which we get $f(a + b + \sqrt{a} f(4b)) = f(b + a + \sqrt{b} f(4a))$. Since $f$ is injective in the case (1), we conclude that $a + b + \sqrt{a} f(4b) = b + a + \sqrt{b} f(4a)$, and hence $\sqrt{a} f(4b) = \sqrt{b} f(4a)$ holds for every pair of non-negative real numbers $a$ and $b$. By substituting $a = 1$ and $b = \frac{x}{4}$, where $x$ is an arbitrary non-negative real number, into the last equation we get $f(x) = \frac{f(4)}{2} \sqrt{x}$. Therefore, we can write $f(x) = k \sqrt{x}$ for some $k > 0$. If we substitute $x = y = 1$ into the functional equation for this choice of $f$, we obtain that $2k = k \sqrt{2} + 2k$, from which we get, since $k > 0$, $2 = \sqrt{2} + 2k$, and therefore, $k = 1$. Thus, we conclude that $f(x) = \sqrt{x}$.

Conversely, if $f(x) = \sqrt{x}$, then we get for any pair of non-negative real numbers $x$ and $y$,
$$\begin{align*} f(x^2) + f(y) &= x + \sqrt{y} \\ f(x^2 + y + x f(4y)) &= \sqrt{x^2 + y + x \sqrt{4y}} = \sqrt{(x + \sqrt{y})^2} = x + \sqrt{y}. \end{align*}$$
So the function $f(x) = \sqrt{x}$ satisfies the given functional equation. Thus, the desired answer in the case (1) is $f(x) = \sqrt{x}$.

Alternate Solution:
The proof of the fact that $f(x)$ is monotone increasing in the wide sense is the same as the preceding solution. We now separate the 2 cases as follows:
(1) When $f(x)$ is injective, i.e., when $f(a) = f(b) \Rightarrow a = b$.
(2) When $f(x)$ is not injective, i.e., when there exist a pair $a, b$ ($a \neq b$) for which $f(a) = f(b)$.
The fact that $f(x) = \sqrt{x}$ for the case (1) can be shown as in the preceding solution.

Let us show that $f(x) = 0$ must hold in the case (2). Choose a pair $a, b$ ($a < b$) for which $f(a) = f(b)$. Let $t$ be an arbitrary non-negative real number and substitute $(x, y) = (t, a)$, $(t, b)$ into the functional equation, then we get
$$\begin{aligned} f(t^2) + f(a) &= f(t^2 + a + t f(4a)), \\ f(t^2) + f(b) &= f(t^2 + b + t f(4b)), \end{aligned}$$
from which we get, in view of the fact that $f(a) = f(b)$,
$$f(t^2 + a + t f(4a)) = f(t^2 + b + t f(4b)).$$
In particular, for $t > 0$ for which $t^2 + a + t f(4a) = b$ (such a $t$ exists since $a < b$) we get $f(b) = f(b')$ if we put $b' = t^2 + b + t f(4b)$. We thus have $f(a) = f(b')$ with $b' = (t^2 + a + t f(4a)) + (b - a) + t(f(4b) - f(4a)) \ge b + (b - a) = 2b - a$ since $f(4b) \ge f(4a)$. Summarizing this we now see that if $f(a) = f(b)$ holds for some pair of non-negative real numbers $a, b$ ($a < b$), then there exists a $b'$ for which
$$f(a) = f(b'), \quad b' \ge 2b - a$$
hold.

Repeating the procedure above to obtain $b'$ from the pair $(a, b)$, we get a sequence $b_1, b_2, \dots$ of positive numbers for which
$$\begin{aligned} \bullet & b_1 = b \\ \bullet & f(a) = f(b_n) \quad (n = 1, 2, \dots), \\ \bullet & b_{n+1} \ge 2b_n - a \quad (n = 1, 2, \dots) \end{aligned}$$
From $b_{n+1} \ge 2b_n - a$ we get $b_{n+1} - a \ge 2(b_n - a)$, from which it follows that $b_n - a \ge 2^{n-1}(b_1 - a)$. Since $b_1 - a > 0$, we see that $b_n$ can be taken to be arbitrarily large by choosing $n$ larger and larger.

Let now $x$ be an arbitrary positive number for which $x \ge a$, then for a sufficiently large $n$ we have $a \le x \le b_n$, and since $f$ is monotone increasing in the wide sense and since $f(a) = f(b_n)$, we get $f(x) = f(a)$. Thus, we have shown that $x \ge a \Rightarrow f(x) = f(a)$. If we take both $x$ and $y$ to be sufficiently large real numbers and substitute these in the functional equation, we get the identity $f(a) + f(a) = f(a)$, from which we conclude that $f(a) = 0$. Consequently, we have $f(x) = f(a) = 0$ if $x \ge a$. Since $f(x)$ is non-negative valued, and is monotone increasing in the wide sense, we also have $0 \le x \le a \Rightarrow 0 \le f(x) \le f(a) = 0$ as well, and therefore, $f(x) = 0$ for all $x \ge 0$.