Maximum Likelihood with Gradient Descent or Coordinate Descent blows up
You correctly computed the gradients, but the optimization step is incorrect. Your gradient is computed at point $(\mu,\sigma)$, and, therefore, it should be used to update $\mu$ and $\sigma$ simultaneously. Correct update rule is
$$
\begin{align}
\mu_{t+1} &amp;\longleftarrow \mu_t + \gamma \left(\frac{\overline{x}}{\sigma^2_t} - \frac{\mu_t}{\sigma^2_t}\right) \\
\sigma_{t+1} &amp;\longleftarrow \sigma_t + \gamma\left(\frac{1}{n\sigma^3_t}\sum^n_{i=1}(x^{(i)} - \mu_{t})^2 - \frac{1}{\sigma_t}\right)
\end{align}
$$
The formula that you have in the question implements the alternating optimization which may or may not converge (it is guaranteed to converge if the gradient norm along the optimization path is at most one).
Therefore, I corrected your code as follows:
# Loop through and update mu and sigma
mus = [mu]
sigmas = [sigma]
for i in range(n_iter):
    # Compute gradients
    mu_grad = (np.mean(x) - mu) / (sigma**2)
    sigma_grad = (np.mean((x - mu)**2)/sigma**3 - 1 / sigma)
    # Update mu and sigma
    mu, sigma = mu + gamma_mu*mu_grad, sigma + gamma_sigma*sigma_grad
    # Store mu and sigma
    mus.append(mu)
    sigmas.append(sigma)

The result converges quickly.