Biquad coefficients using Magnitude (or Phase) Invariance Mapping Method
2nd order MIM LPF results looks great for selected cutoff frequency (fs = 44.1kHz) as seen in plot (frequency range upto 23873Hz):



Coefficients for plotted 2nd order filter:

a [0.437417389844334 0.1220119196463959 -0.04071587087954469]
b [1 -0.6951293456643256 0.2138427842755107]


EDIT:

It looks like there are some problems in magnitude response when fs is increased and/or when fc is reduced. MIM/PIM paper mentions this as an issue.

So, it looks like MIM can be used for to improve 2nd order LPF magnitude response especially when using low fs and high enough fc. Magnitude issues can be bypassed by using higher order filters. Fortunately MIM method works correctly for my purposes.

EDIT: First and the second D(s) given in paper examples are wrong. The correct ones are (5s + 2) and (s + 10) / (s^2 + 4.1s + 0.4). With those all the MIM calculations become right. The third example (PIM) was correct.