BifurcationKit - Error using continuation()

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?

it probably means the jacobian is singular. plase use ‘verbosity=3’ to report where this occurs.

are you sure thisvector field is differentiable?