Simple algorithm for gaussian mixture model

In a gaussian mixture model, I don't really understand how to update the covariance matrix s_j associated to each cluster. Is the following code correct (to see if I've understood how it works):

If I understand, we should repeat the following code algorithm until convergence:

X = {x1, ..., xn} set of data-points
M = {m1, ..., mk} set of k means choosen randomly from X
S = {s1, ..., sk} set of k covariance matrix, each si is choosen to be the identity matrix

for i from 1 to n
{
    for j from 1 to k
    {
        compute Pij = pdf(xi | mj, sj) / SUM_l pfd(xi | ml, sl)
    }
}

for i from 1 to n
{
    for j from 1 to k
    {
        mj = mj + (Pij * xi) / SUM_l Plj // update the mean
        // here how to update sj ? covariance matrix
    }
}