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

