function t = stab221_BW(theta,mu)
t = stab211_BW(theta,mu) + 3*mu*mu*(1 - mu).*exp(-i.*theta).*(exp(i.*theta)-1).*(1 - exp(-i.*theta)).^2;

