function M = linearisedInvader
%The two-patch model linearised for invaders (m with sigma = sigman) around equilibria for residents (n with simga = sigmam)

%output: 2x2 anonymous function array with arguments: (b,del0,del1,eps,lam,mu,n0,n1,sigi,sigr)

syms n0 n1 m0 m1 del0 del1 b lam sigman sigmam
eps = sym('eps');
mu = sym('mu');

f0 = b - del0*(n0+m0) - mu - lam*sigmam;
f1 = b - del1*(n1+m1) - mu - lam*sigmam;

L = eps*sigmam*m0*(sigman*n1 + sigmam*m1)*(f1-f0);

d0 = m0*f0 + (f1 + mu)*m1 - L;
d1 = -m1*mu + L; 

d0d0 = diff(d0,m0); %Finding linearised combonents
d0d1 = diff(d0,m1);
d1d0 = diff(d1,m0);
d1d1 = diff(d1,m1);

M00 = subs(d0d0,{m0,m1},{0,0}); %Setting to (n0r* n1r* 0 0)
M01 = subs(d0d1,{m0,m1},{0,0});
M10 = subs(d1d0,{m0,m1},{0,0});
M11 = subs(d1d1,{m0,m1},{0,0});

M = matlabFunction([M00,M01;M10,M11]);