black scholes - What are the main problems for calculating the implied volatility of in the money American put options?

As stated in the question I have a problem with calculating the implied volatility for in the money put options I have a data set of 2.6 million American style plain-vanilla call and put options. For these options I need to calculate the implied volatility (already wrote a Python script). For call options and out-of-the-money put options, I am able to determine almost 100% of the implied volatilities. However, it does a really bad job for calculating in the money put options (see also the graph below).

For calculating the implied volatility I use the Black-Scholes formulas for plain-vanilla put and call options. My data set only contains American options (instead of the European assumption of BSM) and is downloaded from the Euronext website and freely available (https://www.euronext.com/nl/reports-statistics/derivatives/daily-statistics).

Furthermore, I make use of discrete dividends instead of continuous dividends. Main reasons for this is, that dividends aren't continuous and I need to price barrier options, which are very sensitive for ex-dividend dates. I also discount the expected dividends (ex-dividend) between t and maturity against Eonia rate (same rate as the risk-free for the BSM).

Also this data set contains options from 1-1-2017 until now (negative Eonia rate). Below you see a print-out of the first 10 put options for which I am not able to calculate the implied volatility. The error column is the difference between the quoted price (Quoted_Put) and BSM price (P_BSM). I also added the IV and the quoted price for the call counterpart (same underlying, strike, maturity and quoted date).

My question is pretty simple does some one has an academic or practical suggestion how I can solve this problem? Or maybe someone see a pattern which needs more attention? It surprises me that I do not have a problem with calculating the implied volatility of other options, only in the money put options does give problems.

I checked the formulas twice in Excel and Python (similar results) but maybe did I make a wrong assumption? Or maybe someone encountered a similar problem and can tell me how he or she solved it? If you have any further questions, I am happy to provide extra information, graphs and spreadsheets.

A couple of non-preferable solutions which crossed my mind:

I could use only the out of the money options (in this case I could easily drop most of the errors), for which I need to assume that the implied volatility is the same for the call and the put (well known this is not true and also visible in the spreadsheet). This is not a very strong argument (main argument is that is solves my problem). This solution is also suggested by a contributor on another question (Implied Volatility Calculation for Deep In The Money Calls, Numerical Issues). But it is not really clear why this is rational/academic feasible to do (@Mark Joshi stated: the intrinsic value of the in the money tends to swamp the numerics)

I also could use only call options with the same assumptions as the before mentioned "solution". In an interest rate positive environment an American call and an European call are equal to each other. However, the Eonia is already negative since 2014.

A different risk-free rate (i.e.) Euribor 3M also improves the results but I heard en read that in practice traders use the Eonia rate. So, also not a very strong solution. Furthermore, the Euribor is also negative.