for J=0:Jmax
    hsym=zeros(2*J+1);
    for Kp=-J:J
        m=m+1;
        n=0;
        for K=-J:J
            n=n+1;
            if Kp==K
                hsym(n,m)=F*(J*(J+1)-K^2)+G*K^2;
            elseif Kp==K+2
                hsym(n,m)=H*(0.25*(J*(J+1)-K*(K+1))*(J*(J+1)-(K+1)*(K+2)))^(0.5);
            elseif Kp==K-2
                hsym(n,m)=H*(0.25*(J*(J+1)-K*(K-1))*(J*(J+1)-(K-1)*(K-2)))^(0.5);
            end
        end
    end
    [evrp,evlp]=eig(hsym);
    %%evlp=evlp';
    for e=1:(2*J+1)
        evl(J^2+e,J^2+e)=0.5*(A+C)*J*(J+1)+0.5*(A-C)*evlp(e,e);
        for ee=1:(2*J+1)
            evr(J^2+e,ee)=evrp(e,ee);
        end
    end
end