Here is code that adjusts r at the frequencies near each pole pair, until -0.9<=H(w)<=-1.0. Since there are three pole pairs, there are three values for r: one for each pole pair. The plot below shows Hinit (=H before adjusting the three r's), and final H, after adjusting three r's. The later r's are affected by the r's already adjusted. If you change the order of adjusting the r's (3,2,1 versus 1,2,3), you get slightly different values for the three r's, but the transfer function is not noticeably different. Code to make the figure is attached.
In your original code, you have
H = 20*log10(r*abs(Freq));
H = 20*log10(r*abs(freqz(b,a,phi)));
Why do you multiply abs(Freq) by r, and abs(freqz()) by r? That seems unjustified - in the sense that, if you filter a signal with
then the transfer function from x to y will not inclue the factor r that you are including in H. Also, since there are three r's (one for ech pole pair), it is not obvious which r shoudl be used. Therefore I have removed the factor of r from Hinit and from H.