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