clear all

Omega=linspace(0.95,2,801);
EP=linspace(0,0.9,501);

for ie=1:length(EP);
    e=EP(ie);
    for io=1:length(Omega);
        om=Omega(io);
        o1=om+e;
        o2=om-e;
        c1=cos(pi*o1);s1=sin(pi*o1);
        c2=cos(pi*o2);s2=sin(pi*o2);
        PHI1=[c1, s1/o1;-o1*s1,c1];
        PHI2=[c2, s2/o2;-o2*s2,c2];
        PHI=PHI1*PHI2;
        sigma=max(abs(eig(PHI)));
        if sigma<1.0000001;
            lambda(io,ie)=-10;
        else;
        lambda(io,ie)=log(sigma)/(2*pi);
        end; 
    end;
end;
    v=linspace(0.0001,0.2,11);
    figure(32)
    contourf(Omega,EP,lambda',v);colormap jet;colorbar
    xlabel('\omega');ylabel('\epsilon')