Generally this is accomplished (i.e. simplified) using matrix notation. Let $X$ be a matrix whose rows are your data such that the columns are your dimensions. Further, let $y$ be a vector of your response variable and $\beta$ be a vector of the regression coefficients (i.e. what you are solving for.
Append a column of 1's to the left hand side of the matrix $X$ to construct what is called the "design" matrix. This allows you to represent your data as:
$y = X \beta$
Then the normal equations are:
$X^T y = X^TX \beta$
And so the solution for the coefficients is:
$ \hat{\beta}_{OLS}= (X^TX)^{-1} X^T y$
Now, although this analytic solution exists, it's often numerically unstable. The size of the matrix $X^TX$ scales with the square of the number of records in your data, and so inverting it is problematic (and the matrix is often "ill-conditioned," which I won't get into but basically means inverting it numerically is hard).
All this is to say: this is one simple way to solve multiple regressions, but you're probably much better off using a software solution maintained by people who care a lot about doing this well and efficiently. I found this library, which looks like it has regression capabilities, but I think it's just using the process I described above in which case it might not be able to accomplish regressions over a lot of data. But you can still give it a shot and see if it works for your use case.
If you really want to roll the solution yourself and you find the matrix inversions to be problematic, this article describes matrix decompositions that will help you get around some of the issues I described.