$y_i = \beta_0 + \beta_1x_{i1} + ... + \beta_px_{ip} + \epsilon_i, i = 1,2,\ldots,n$
where $ \epsilon_i \sim N(0,\sigma^2)$ are iid.
In this case, I am trying to find the MLE of $\boldsymbol{\beta}$ and $\sigma^2$, where $\boldsymbol{\beta} = (\beta_0,\beta_1,...,\beta_p)^\top$.
I am aware that to find MLE of a simple linear regression, I could find the Log-likelihood, before differentiating it by $\frac{dL}{d\beta_0}$ and $\frac{dL}{d\beta_1}$ and equating them to $0$ to find $\hat{\beta_0}, \hat{\beta_1}$ and $\hat{\sigma^2}$.
By extending the method for finding MLE of simple linear regression, I tried the following method:
$\ell(\boldsymbol{\beta},\sigma^2)= \log(P(y|\boldsymbol{\beta},\sigma^2)) =-\frac{n}{2}\log(2\pi)-\frac{n}{2}\log(\sigma^2)-\frac{1}{\sigma^2}\sum^n_{i=1}(y_i-\beta_0-\beta_1x_{i1}- \ldots -\beta_px_{ip})^2$
However, I am stuck as when I try to differentiate the log-likelihood w.r.t. $\beta_0, \beta_1, \ldots, \beta_p$.
Am I on the right track? How am I supposed to proceed?