Suppose your tree is splitting some continuous variable $x.$ In your data, you have $M$ observations, each with a specific value for this particular column, that already falls into the given branch before splitting, i.e. $\{x_1, x_2, \dots, x_M\}$. Note, that $M$ isn't necessarily the original size of the data, because this may be one of many splits that have already occurred, and not all the data necessarily falls into the branch where this new split is taking place.
For the $M$ points, after the split they will take on values $\{ \hat{y_1}, \dots, \hat{y_M} \}$, each of which will be one of two possible values: $w$ or $z$. We'll say that $\hat{y_i} = w$ if $x_i < \lambda$, and $\hat{y}_i = z$ otherwise. Thus $x = \lambda$ is the splitting point.
Suppose you are trying to minimize some loss function, $L(\textbf{y}, \hat{\textbf{y}}) = \sum_i l(y_i, \hat{y_i}).$ This may be the sum of squared error, the cross-entropy, or something else depending on whether you're doing regression, classification, etc. For example, if your loss function is the sum of square error (assuming you're doing regression), the post-branch value of your loss function for all $M$ values in the current branch is:
$$
L = \sum_{x_i < \lambda} \left ( y_i - w \right)^2 + \sum_{x_i \ge \lambda} \left(y_i - z \right)^2.
$$
The goal is to minimize this with repsect to $\lambda$, $w$, and $z$ simultaneously.
It's not easy to see how to minimize this right off the bat, because it's a discontinuous function with respect to $\lambda$ (changing it changes which $x$'s are in each region). For a given $w$ and $z$, $L$ is a piecewise flat function of $\lambda$. That is, if there are no values in the observation set between $x_k$ and $x_l$, then $L$ is a constant in $[x_k, x_l).$ Thus, for convenience, $\lambda$ can be taken to be any one of points in the observed values $\{x_1, \dots, x_M \},$ because any slight shift in $\lambda$ that doesn't change the region distribution of the $x$'s will not change $L$. It's easy to see that, for a given value of $\lambda$, the optimal values of $w$ and $z$ are the means of the $y$ values in each respective region. Thus, to optimize $w, z, \lambda$, we can simply run $\lambda$ through all values in $\{x_1, \dots, x_M\}$, assigning $w$ and $z$ to the appropriate means in each case, and evaluate $L$. We choose the values that minimizes $L.$
This tells you how to find a branching point for a given feature $x$. However, at each branch, you want to minimize not only the branching point and values of each branch, but the variable that you want to branch. Since the minimization for each variable is a quick process, it suffices to do the minimization for all variables (or all of them in the subset being considered, if you're subsampling the columns to avoid overfitting), then pick the one that gives the best minimization (biggest decrease in loss function).