Pre-Warping is not working on bilinear transform New here and I just wanna ask the question of what happens when your pre warping doesnt work on your bilinear transform?
I am assuming my isnt working as part of my research I have read that the higher the frequency is the more non linear it becomes, so I am assuming that my frequency I chose just has to many non linearities that it wouldn't ever get re-aligned to the continuous function critical frequency. It tries  but doesnt get there.
Sampling Frequency = 44.410kHz
Nyquist Frequency = 22.205kHz
Critical Frequency = 20kHz
Ratio = Fs/Nyquist frequency =  20/22.205 = 90.1%
I noticed the closer you are critical frequency is to your Nyquist frequency when applying the bilinear transform it gets ugly. And it looks like the pre warping can only do so much.
The only way I see around this if I just changed my Sampling frequency to like 96kHz which would make my Nyquist frequency = 48kHz would be then the ratio = 42% this would ease the constraint of the non-linearities and have the pre warping affect work much better, but this raises the question is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz? Because the way this is going I dont see it.

CODE:
s = tf('s');

Fs = 96e3;
Ts = 1/Fs;
Fc = 20e3;

GP = -3;
GS = -20;
WP = 2*pi*Fc;
WS = 251327.412;

n = ceil(log10((10^(-GS/10)-1)/(10^(-GP/10)-1))/(2*log10(WS/WP)));

WC_1 = WP/((10^(-GP/10)-1)^(1/(2*n))); 
WC_2 = WS/((10^(-GS/10)-1)^(1/(2*n))); 

LP = 1/((s/WC_1)^4+2.613*(s/WC_1)^3+3.414*(s/WC_1)^2+2.613*(s/WC_1)+1);
LP_2 = 1/((s/WC_2)^4+2.613*(s/WC_2)^3+3.414*(s/WC_2)^2+2.613*(s/WC_1)+1);

LP_44_410kHz = c2d(LP,1/44.410e3,['Method','tustin','PrewarpFrequency',WP]);
LP_96kHz = c2d(LP,Ts,['Method','tustin','PrewarpFrequency',WP]);

options = bodeoptions;
options.FreqUnits = 'Hz';
options.Title.String = 'With Pre-Warping';
bode(LP_44_410kHz,LP,LP_96kHz,options);
is it possible to have a digital filter with a sampling frequency of 44.410kHz and a Fc = 20kHz?

Yes it is possible... .
I use MIM (Magnitude Invariance Method) for this. Here's plot showing 2nd order LPF responses for fc=20kHz-25kHz (1kHz step) @ 44.kHz sampling:

MIM method works well for even lower sample rate implementations but, there are also issues not sorted out yet.
Here's Matlab implementation of the MIM and PIM (Phase Invariance Method) method.