I would just like to add a practical example to Reto Koradi's answer.
Let's assume you already have a 4x4 model matrix and want to use it to transform normals as well. You can start by deducing scale in each axis by taking length of the 3 first columns of that matrix. If you now divide each column by its corresponding scaling factor, the matrix will no longer affect model's scale, because the basis vectors will have unit length.
As you pointed out, normals have to be scaled by the inverse of the scale in each axis. Fortunately, we have already derived the scale in the first step, so we can divide the columns again.
All that effectively means that if you want to derive transform matrix for normals from your model matrix, all you need to do is to divide each of its first three columns by their lengths squared (which can be rewritten as dot products). In GLSL you would write:
mat3 mat_n = mat3(mat_model);
mat_n[0] /= dot(mat_n[0], mat_n[0]);
mat_n[1] /= dot(mat_n[1], mat_n[1]);
mat_n[2] /= dot(mat_n[2], mat_n[2]);
vec3 new_normal = normalize(mat_n * normal);