Hello, I am new to Julia and trying to use BifurcationKit to do some continuation analysis. However, I have yet to get it working for my system of ODEs. Having debugged a few errors, I am now getting this one that I cannot work out how to resolve:
LinearAlgebra.SingularException(9)
Here is my code:
using BifurcationKit
function model(X,p)
(;b1, f, d, alpha, gamma, rho, beta, LR, LT) = p
P0, PS, PI, P00, PS0, PI0, PSI, PSS, PII = X
b2 = f*b1
if P0>0
q_S0 = PS0/P0
q_I0 = PI0/P0
else
q_S0 = 0
q_I0 = 0
end
if PS>0
q_IS = PSI/PS
else
q_IS = 0
end
if PI>0
q_SI = PSI/PI
else
q_SI=0
end
dX1 = LR*(-(b1*q_S0+b2*q_I0)*P0)-(1-LR)*P0*(b1*PS+b2*PI)+(d+alpha)*PI+d*PS
dX2 = LR*(b1*q_S0+b2*(1-rho)*q_I0)*P0+(1-LR)*P0*(b1*PS+b2*(1-rho)*PI)-LT*beta*q_SI*PI-(1-LT)*beta*PS*PI-d*PS+gamma*PI
dX3 = LR*b2*rho*q_I0*P0+(1-LR)*b2*rho*P0*PI+LT*beta*q_SI*PI+(1-LT)*beta*PS*PI-(d+alpha+gamma)*PI
dX4 = 2*(d+alpha)*PI0+2*d*PS0-2*LR*(3/4)*(b1*q_S0+b2*q_I0)*P00-2*(1-LR)*(b1*PS+b2*PI)*P00
dX5 = LR*(3/4)*(b1*q_S0+b2*(1-rho)*q_I0)*P00+(1-LR)*(b1*PS+b2*(1-rho)*PI)*P00+gamma*PI0+d*PSS+(d+alpha)*PSI-d*PS0-LT*beta*(3/4)*q_IS*PS0-(1-LT)*beta*PI*PS0-LR*(b1/4+(3/4)*(b1*q_S0+b2*q_I0))*PS0-(1-LR)*(b1*PS+b2*PI)*PS0
dX6 = LR*b2*(3/4)*rho*q_I0*P00+(1-LR)*b2*rho*PI*P00+LT*(3/4)*q_IS*beta*PS0+(1-LT)*beta*PI*PS0+d*PSI+(d+alpha)*PII-(d+alpha+gamma)*PI0-LR*(b2/4+b2*(3/4)*q_I0+b1*(3/4)*q_S0)*PI0-(1-LR)*(b1*PS+b2*PI)*PI0
dX7 = LR*(b1*(3/4)*q_S0+b2*(1-rho)*(1/4+(3/4)*q_I0))*PI0+(1-LR)*(b2*(1-rho)*PI+b1*PS)*PI0+LR*(3/4)*b2*rho*q_I0*PS0+(1-LR)*b2*rho*PI*PS0+LT*(3/4)*beta*q_IS*PSS+(1-LT)*beta*PI*PSS+gamma*PII-(2*d+alpha+gamma)*PSI-LT*beta*(1/4+(3/4)*q_IS)*PSI-(1-LT)*beta*PI*PSI
dX8 = 2*LR*(b1/4+b1*(3/4)*q_S0+b2*(3/4)*(1-rho)*q_I0)*PS0+2*(1-LR)*(b1*PS+b2*(1-rho)*PI)*PS0+2*gamma*PSI-2*d*PSS-2*LT*(3/4)*beta*q_IS*PSS-2*(1-LT)*beta*PI*PSS
dX9 = 2*LR*b2*rho*(1/4+(3/4)*q_I0)*PI0+2*(1-LR)*b2*rho*PI*PI0+2*LT*beta*(1/4+(3/4)*q_IS)*PSI+2*(1-LT)*beta*PI*PSI-2*(d+alpha+gamma)*PII
return [dX1,dX2,dX3,dX4,dX5,dX6,dX7,dX8,dX9]
end
parameters = (b1 = 0.5, f = 0.5, d=1., alpha=0.5, gamma=1.5, rho=0.75, beta=10., LR=0.5, LT=0.5)
x0 = [0.42372359, 0.45432498, 0.12195143, 0.1929949 , 0.18224381,0.04848488, 0.04775254, 0.22432863, 0.02571401]
prob = BifurcationProblem(model, x0, parameters, (@optic _.b1))
br = continuation(prob, PALC(), ContinuationPar(p_min = 0., p_max = 10.), bothside = true, verbosity = 3)
Does anyone have any thoughts on what the problem might be?