%Code for section 4.3.2

a = f(1:19800);%init data representing the population of species A
b = r(1:19800);%init data representing the population of species B

%%FOXES matrix equation coefs
m11 = -sum(a*a')
m12 = sum(a.*a.*b)
m21 = -m12
m22 = sum(a.*a.*b.*b)

z1 = (f(2:19801)-a)*a'
z2 = sum((f(2:19801)-a).*a.*b)

M = [m11 m12; m21 m22];%now compute the matrix by assembling the components
R1 = [z1;z2]%this is the RHS of equation for (g_1, c_1)

cg1 = inv(M)*R1%solve for (g_1, c_1)


%%rabbits equation coefs
m11 = sum(b*b')
m12 = -sum(b.*b.*a)
m21 = -m12
m22 = -sum(a.*a.*b.*b)

z1 = (r(2:19801)-b)*b'
z2 = sum((r(2:19801)-b).*a.*b)

M = [m11 m12; m21 m22];
R1 = [z1;z2]

cg2 = inv(M)*R1%solve for (g_2, c_2)