clear variables
syms omega
thetas = [
    pi/6;
    pi/4;
    pi/2;
    2*pi/3;
    5*pi/6
];

Zin = exp(1j*omega);
for iter = 1:length(thetas)
    BZ = 1/(1 - 2*cos(thetas(iter))*(Zin^(-1)) + (Zin^(-2)));
    dBZ = diff(BZ);
    Zval = solve(dBZ == 0, omega);
    Zcalc = exp(1j*Zval);
    BZOut(iter) = double(1/(1 - 2*cos(thetas(iter))*(Zcalc^(-1)) + (Zcalc^(-2))));
end