Bilinear Transformation Comparison
i really don't wanna slog through your MATLAB code.  is this what you're trying to prove? :

let $H_\text{a}(s)$ be some analog (or "continuous-time") LTI transfer function on the s-plane and $H_\text{d}(z)$ be some digital (or "discrete-time") LTI transfer function in the z-plane.

$$ H_\text{d}(z) = H_\text{a}(s) \Bigg|_{s=\frac{2}{T}\frac{z-1}{z+1}} $$

do you want to prove that

$$ H_\text{d}(e^{j\omega}) = H_\text{a}(j\Omega) \Bigg|_{\Omega=\frac{2}{T}\tan\left(\frac{\omega}{2}\right)} $$

is that what you want proven?