In multivariate linear regression (linear regression with more than one variable) the model is $yi = b_0 + b_1x_{1i} + b_2x_{2i} + ...$ , and so on. But how is the $b_n$ value calculated iteratively? Can it be calculated non-iteratively? What is the intuition behind using that method to calculate $b_2$?
It is calculated the same way $b_1$ is calculated.
Nearly following your notation, say your multiple linear regression function is
$H(X_i) = b_0 + b_1x_{1,i} + ...+ b_nx_{n, i}$
for data instance $X_i=x_{1,i},...,x_{n, i}$ and weights $b_0,...,b_n$.
And say your error function is $E(X,Y) = \sum_i(H(X_i)-Y_i)^2$
where $X$ is the collection of all data points $X_i, Y_i$.
From your error function $E$, for whatever weights you have (with a gradient based method), calculate the partial derivative $\partial E /\partial b_i$ and use this to update all of your weights at once in each iteration of your optimization routine.
put your b_1, ... ,b_n coefficients into a vector b
put all your x_{ij} into a matrix X
then all components of b are calculated at the same time with this equation
b = (H(X) - Xb)^T (H(x) - Xb)
But this calculation (estimation) is only consistent (search consistent estimator) when certain assumptions are present (read here https://en.wikipedia.org/wiki/Ordinary_least_squares).