I'm trying to model some data using the following distribution: \begin{align} r &\sim \mathcal{N}(\mu, \sigma^2) \nonumber\\ \mu &= m_0 + m_1 d \nonumber\\ \sigma^2 &= s_0+s_1 d\nonumber\\ p(r_1, r_2, \dots, r_N) &= \prod_{i=1}^N \frac{1}{\sqrt{2\pi\sigma^2_i}}e^{\frac{1}{-2\sigma^2_i}(r_i-\mu_i)^2} \nonumber\\ \mathcal{L} &= -\frac{1}{2}\log\big(p(r_1, r_2, \dots, r_N)\big)\nonumber\\ &= -\frac{1}{2}\sum_{i=1}^N\log\Big(\frac{1}{\sqrt{2\pi\sigma^2_i}}e^{\frac{1}{-2\sigma^2_i}(r_i-\mu_i)^2}\Big) \nonumber\\ &= \sum_{i=1}^N\Big[\log\big(2\pi\sigma^2_i\big) + \frac{(r_i-\mu_i)^2}{\sigma^2_i}\Big]\nonumber \end{align}
I was wondering if it is possible to maximise this likelihood analytically to obtain fits for $m_0, m_1, s_0, s_1$. I know that there exists a solution if $\sigma^2$ is constant (simple OLS) but what if it isn't? (I tried the first few steps below but am stuck on the last line)
\begin{align} \frac{\partial\mathcal{L}}{\partial s_0} &= \sum_{i=1}^N \frac{\partial}{\partial \sigma^2_i}\frac{\partial\sigma^2_i}{\partial s_0}\Big[\log\big(2\pi\sigma^2_i\big) + \frac{(r_i-\mu_i)^2}{\sigma^2_i}\Big]\nonumber\\ &=\sum_{i=1}^N\frac{1}{\sigma_i^2}-\frac{(r_i-\mu_i)^2}{\sigma_i^4}=0 \nonumber\\ \end{align}