title: Error in lsoda(y, times, func, parms, ...) : illegal input detected before taking any integration steps
body: <p>I am trying to estimate and optimize parameters with time dependent transmission rates in the infectious disease dynamics model. Since the incidence data is periodic, assume betaH=aH*(1+bH<em>cos(2</em>pi*t/12)) according to the literature. But under real data simulation, I get the error: illegal input detected before taking any integration steps.
  I'm not very familiar with differential equation solving, so please tell me, how it should be modified.</p>
  <pre><code>library(deSolve)
  
  # SEIARWF differential equation
  SEIARWF &lt;- function(t, state, parameters) {
    with(as.list(c(state, parameters)), {
      dS=0.001048*N+0.00595*R-(aH*(1+bH*cos(2*pi*t/12)))*S*(I+A)/N-
                              (aW*(1+bW*cos(2*pi*t/12)))*S*W/N-
                              (aF*(1+bF*cos(2*pi*t/12)))*S*FD/N-0.000582*S
      dE=(aH*(1+bH*cos(2*pi*t/12)))*S*(I+A)/N+
         (aW*(1+bW*cos(2*pi*t/12)))*S*W/N+
         (aF*(1+bF*cos(2*pi*t/12)))*S*FD/N-
         0.75*E-0.000582*E
      dI=0.1*0.75*E-1*I-0.000582*I
      dA=0.9*0.75*E-1*A-0.000582*A
      dR=1*I+1*A-0.00595*R-0.000582*R
      dW=3014000000000*I+1507000000000*A-1.3333*W
      dFD=3014000000000*I+1507000000000*A-1.3333*FD
      N=S+E+I+A+R
      
      list(c(dS, dE, dI, dA, dR, dW, dFD, N)) 
    })
  }
  
  # state
  N0 &lt;- 1292270000                   
  Infected &lt;- c(1429,1813,2227,1752,1541,1467,1249,1188,1022,973,1052,731,1539,1604,
                2190,1838,1327,1183,1038,1052,942,997,1045,786,1704,2219,2768,2116,
                1818,1490,1228,1207,1042,1179,1213,1023,2183,2039,2616,2176,1842,1643,
                1571,1612,1334,1306,1292,963,1843,1871,2726,2112,1786,1425,1305,1226,
                1061,1096,1097,977,1930,2444,2730,2048,1735,1487,1338,1424,1157,1195,
                1237,1550,2288,2304,3076,2567,2180,1719,1584,1461,1421,1468,1684,1930,
                2342,2930,4178,3270,2459,2118,2010,2116,1871,1872,2062,1974,2338,3180,
                3557,2662,2354,1889,1843,1906,1737,1763,1959,2083,3076,3145,3628,2772,
                2241,1794,1874,1899,1723,1764,1960,2026,2655,3018,3140,2510,2052,1836,
                1958,2054,1908,1849,1883,2125,2893,2301,3032,2668,2165,1956,1991,1931,
                1792,1982,2131,2327,2497,2468,3318,2676,2410,2143,2094,2105,1966,1885,
                2065,2295,2306,2752,3141,2568,2559,2397,2390,2439,2049,2006,2191,2216,
                2811,2441,3163,2651,2489,2160,2324,2249,1928,1863,2304,2220)
  I0 &lt;- Infected[1]                  
  A0 &lt;- Infected[1]*9                
  E0 &lt;- 19054                        
  R0 &lt;- 303166542                   
  S0 = N0 - E0 - I0 - A0 - R0        
  W0 &lt;- 5169000000000000             
  FD0 &lt;- 11080000000000000           
  initials &lt;- c(S = S0, E = E0, I = I0, A =A0, R = R0, N = N0, W = W0, FD = FD0)  
  
  # time
  t &lt;- 0:(length(Infected)-1)
  
  # RSS optim
  RSS &lt;- function(parameters) { 
    names(parameters) &lt;- c(&quot;aH&quot;, &quot;bH&quot;, &quot;aW&quot;, &quot;bW&quot;, &quot;aF&quot;, &quot;bF&quot;)
    out &lt;- ode(y = initials, times = t, func = SEIARWF, parms = parameters)
    fit &lt;- out[ , 4]
    sum((Infected - fit)^2)
  }
  
  Opt &lt;- optim(par=c(0.0005, 0.5, 0.0102, 0.5, 0.0429, 0.5), 
               fn = RSS, method = &quot;L-BFGS-B&quot;,
               lower = c(0, -1, 0, -1, 0, -1), 
               upper = c(Inf, 1, Inf, 1, Inf, 1),
               control = list(), hessian = TRUE)
  Opt
  Opt_par &lt;- setNames(Opt$par, c(&quot;aH&quot;, &quot;bH&quot;, &quot;aW&quot;, &quot;bW&quot;, &quot;aF&quot;, &quot;bF&quot;))
  Opt_par
  </code></pre>