syms x y a b c d
vars = [x, y]
eqs = [x*(1-x-1/(x+d)*a*y), y*(-c+1/(x+d)*b*x)];
[xc, yc] = solve(eqs(1), eqs(2));
[xc, yc]
A = jacobian(eqs, vars);
disp('Matrix of linearized system:')
subs(A, vars, [xc(1), yc(1)])
disp('eigenvalues:')
eig(ans)
disp('Matrix of linearized system:')
subs(A, vars, [xc(2), yc(2)])
disp('eigenvalues:')
eig(ans)
disp('Matrix of linearized system:')
subs(A, vars, [xc(3), yc(3)])
disp('eigenvalues:')
eig(ans)