%% Gives back 2*2 stability matrix
function [SM]=Stability_Matrix(d,gs,S,a1,a2,lambda1,lambda2)
epsilon=6-d;

%% Common factors
I2=0;
I2_dg=zeros(1,2);
for X1=1:2
    for X2=1:2
        I2=I2+S(X1,X2)*gs(X1)*gs(X2);
        I2_dg(X1)=I2_dg(X1)+2*S(X1,X2)*gs(X2);
    end
end

I4=0;
I4_dg=zeros(1,2);
for X1=1:2
    for X2=1:2
        for X3=1:2
            for X4=1:2
                for X5=1:2
                    I4=I4+S(X1,X5)*a1(X5,X2,X3,X4)*gs(X1)*gs(X2)*gs(X3)*gs(X4);
                    I4_dg(X1)=I4_dg(X1)+S(X1,X5)*a1(X5,X2,X3,X4)*gs(X2)*gs(X3)*gs(X4);
                    I4_dg(X1)=I4_dg(X1)+3*S(X2,X5)*a1(X5,X1,X3,X4)*gs(X2)*gs(X3)*gs(X4);
                end
            end
        end
    end
end

I3=zeros(2,1);
I3_dg=zeros(2,2);
for X=1:2
    for X1=1:2
        for X2=1:2
            for X3=1:2
                I3(X)=I3(X)+a1(X,X1,X2,X3)*gs(X1)*gs(X2)*gs(X3);
                I3_dg(X,X1)=I3_dg(X,X1)+3*a1(X,X1,X2,X3)*gs(X2)*gs(X3);
            end
        end
    end
end

I5A=zeros(2,1);
I5A_dg=zeros(2,2);
for X=1:2
    for X1=1:2
        for X2=1:2
            for X3=1:2
                for X4=1:2
                    for X5=1:2 
                        I5A(X)=I5A(X)+a2(X,X1,X2,X3,X4,X5)*gs(X1)*gs(X2)*gs(X3)*gs(X4)*gs(X5);
                        I5A_dg(X,X1)=I5A_dg(X,X1)+3*a2(X,X1,X2,X3,X4,X5)*gs(X2)*gs(X3)*gs(X4)*gs(X5);
                        I5A_dg(X,X1)=I5A_dg(X,X1)+2*a2(X,X2,X3,X4,X1,X5)*gs(X2)*gs(X3)*gs(X4)*gs(X5);
                    end
                end
            end
        end
    end
end

I5B=zeros(2,1);
I5B_dg=zeros(2,2);
for X=1:2
    for X1=1:2
        for X2=1:2
            for X3=1:2
                I5B(X)=I5B(X)+a1(X,X1,X2,X3)*gs(X1)*gs(X2)*I3(X3);
                I5B_dg(X,X1)=I5B_dg(X,X1)+2*a1(X,X1,X2,X3)*gs(X2)*I3(X3);
                I5B_dg(X,:)=I5B_dg(X,:)+a1(X,X1,X2,X3)*gs(X1)*gs(X2)*I3_dg(X3,:);
            end
        end
    end
end

%% Derivatives of anomalous dimension and beta functions
self_prop=(-epsilon/2);
self_prop=self_prop+((1/4)+epsilon*lambda1)*I2;
self_prop=self_prop+(-(11/144))*(I2^2);
self_prop=self_prop+((1/6)+2*lambda1+(1/2)*lambda2)*I4;
self_prop_dg=zeros(1,2);
self_prop_dg=self_prop_dg+((1/4)+epsilon*lambda1)*I2_dg;
self_prop_dg=self_prop_dg+(-(11/72))*(I2*I2_dg);
self_prop_dg=self_prop_dg+((1/6)+2*lambda1+(1/2)*lambda2)*I4_dg;

% Beta functions
SM=self_prop*eye(2)+gs*self_prop_dg;
SM=SM+(-1+epsilon*lambda2)*I3_dg;
SM=SM+((7/24)-2*lambda1-(1/2)*lambda2)*(I3*I2_dg+I2*I3_dg);
SM=SM-(1/2)*I5A_dg;
SM=SM+(-(3/4))*I5B_dg;

end