function [f,eqs,stability] = pairwiseInvasability(b,mu,del0,del1,lam,eps,Sigma)
%Finds invasion fitness for all pairs Sigma*Sigma for the two-patch model in (up to) two resident steady states
%and the adaptive equilibria and their stability

%input:
%b,mu,del0,del1,lam,eps: model parameters
%Sigma: vector of sigmas used for both resident values and invader values.


M = linearisedInvader; %loads linearised invader fitness matrix function
numb = length(Sigma);
f = zeros(numb,numb,2);

N0 = zeros(numb,2);
N1 = zeros(numb,2);
for r = 1:numb %loops for every value of resident sigma
    [eq0,eq1,stability] = equilibria(b,mu,del0,del1,lam,eps,Sigma(r)); %finds eq. of resident-only system
    nStab = sum(stability == -1);
    N0(r,1:nStab) = eq0(stability == -1); %set initial state to equilibria
    N1(r,1:nStab) = eq1(stability == -1);

    [N1(r,1:nStab),ord] = sort(N1(r,1:nStab)); %lowest equilibrium as index 1
    N0(r,1:nStab) = N0(r,ord);
    
    for j = 1:2
        for i = 1:numb
            [~,lambda] = eig(M(b,del0,del1,eps,lam,mu,N0(r,j),N1(r,j),Sigma(i),Sigma(r))); %invasion fitness is (log) max. eigenvalue of invader matrix 
            f(i,r,j) = max(diag(lambda));
        end
        if nStab == 1
            f(:,r,2) = f(:,r,1); %if only one eq. then both indeces share value
            break;
        end
    end

    disp([num2str(r) ' / ' num2str(numb)])
end

eqs = [];
stability = [];
for j = 1:size(N0,2)
    dfd = zeros(numb-1,1);
    for i = 2:numb-1
        [~,lu] = eig(M(b,del0,del1,eps,lam,mu,N0(i,j),N1(i,j),Sigma(i+1),Sigma(i))); %Find fitness for invaders with *just* larger sigma
        [~,ld] = eig(M(b,del0,del1,eps,lam,mu,N0(i,j),N1(i,j),Sigma(i-1),Sigma(i))); %Find fitness for invaders with *just* smaller sigma
        dfd(i) = max(diag(ld))-max(diag(lu)); %Difference in invasion fitness above and below the diagonal
    end
    dfd(1) = [];
    fisx = dfd(1:end-1).*dfd(2:end);
    eqidx = find(fisx<0); %set equilibria to where dfd changes sign
    eqs = [eqs Sigma(eqidx + 1)]; %dynamically adds eq. to output
    ddfd = diff(dfd);
    stab = ddfd(eqidx);
    stability = [stability;stab + 1 - 2*(stab>0)]; %find stability from second order (numerical) derivative
end