Delay Differential Equation Interruptions

I’m trying to simulate a Covid model using a delay differential equation in DifferentialEquations.jl but I’m having a bit of difficultly getting it to work. I’m not sure If I’m applying the methods correctly.

The model is as follows.

function covid_model_delay(beta = 0.2,
                           E_E = 0.3,
                           E_Q = 0.2,
                           E_J = 0.36,
                           mu = 0.000034,
                           z = 0.06,
                           gamma_1 = 0.4,
                           gamma_2 = 0.5,
                           kappa_1 = 0.1,
                           kappa_2 = 0.125,
                           sigma_1 = 0.0337,
                           sigma_2 = 0.0386,
                           d_1 = 0.0079,
                           d_2 = 0.0068)


    S_0 = 4000000.0         # Susceptible
    E_0 = 6.0               # Asymptomatic
    Q_0 = 0.0               # Quarantined 
    I_0 = 1.0               # Symptomatic
    J_0 = 0.0               # Isolated
    R_0 = 0.0               # Recovered
    N_0 = sum(SA[S_0, E_0, Q_0, I_0, J_0, R_0])   # Population

    h(p,t) = ones(7)
    tau = 1
    lags = [tau]

    function F(u, h, p, t)
        histI = h(p, t-tau)[4]

        S, E, Q, I, J, R, N = u
        
        dS = 136.0 - (S*(beta*histI + E_E*beta*E + E_Q*beta*Q +E_J*beta*J)) / N - mu*S
        dE = z + (S*(beta*histI + E_E*beta*E + E_Q*beta*Q +E_J*beta*J)) / N - (gamma_1 + kappa_1 + mu)*E
        dQ = gamma_1*E - (kappa_2 + mu)*Q
        dI = kappa_1*E - (gamma_2 + sigma_1 + d_1 mu)*histI
        dJ = gamma_2*histI + kappa_2*Q - (sigma_2 + d_2 + mu)*J
        dR = sigma_1*histI + sigma_2*J - mu*R
        dN = dS + dE + dQ + dI + dJ + dR

        return SA[dS, dE, dQ, dI, dJ, dR]
    end 

    u0 = SA[S_0, E_0, Q_0, I_0, J_0, R_0]
    tspan = (0,700)
    prob = DDEproblem(F, u0, h, tspan; constand_lags=lags)
    alg = MethodOfSteps(Rodas4(), constrained=true)
    solution = solve(prob, alg)
    return solution 
end 

The function runs and generates a sensible solution for tau=1 as specified, but I can’t run the model for tau greater than 1. Ideally, I want a time lag of 20. The interruption states that a larger max_iter may be required but that doesn’t seem to help. I think I’m made a mistake somewhere along the line.

Any advice would be greatly appreciated.

Is the model divergent for larger tau?