I have also tried all the options of the c2d function with no better results. Converting from transfer function to state space before doing the continuous to digital conversion and convert it back to transfer function gives the same result.
I have also tried to convert it to state space, convert it to digital, convert it to frd, and then identify it in digital --> same result.
Maybe a solution would be to decompose this filter into several filters of smaller order, convert it, and multiply them in the digital domain. I know that c2d is working fine for small order models, like biquadratic filters, any idea on how to do this ?