Creating IIR filter from absorption data

Hello all. I'm working on a FDTD simulation, and I need to create IIR filters to represent boundary conditions.

What I want as input is frequency-dependent absorption coefficients, alpha, between 0 and 1, and reflection factors are calculated as sqrt(1-alpha). Then I want a filter R(z) which obeys R(z)<= 1.

I've tried using invfreqz, but this either returns unstable filters or unsuitable filters (e.g. R(z) > 1).

Example:

f = [63 125 250 500];
alpha = [0.1 0.2 0.25 0.5];
N = 3; % Filter order
w = f/fs*2*pi; % Frequency in radians
R = sqrt(1-alpha);
% Find R(z) coefficients:
[b,a] = invfreqz(R,w,N,N,[],100,1); 
% or
[b,a] = invfreqz(R,w,N,N);

Only with certain combinations I get usable results.

Any suggestions?

Regards, Jakob V