title: array must not contain infs or NaNs when changing np.linspace
body: <p>I am performing sensitivity analysis in Python and want to produce scatterplots to show my system's results. I have been getting results with <code>np.linspace(0,100,100)</code>, but I want to increase my time to be <code>np.linspace(0,1000,100)</code> to see more results in the scatterplot. For some unknown reason, when I increase my time interval, my system does not evaluate (even though time should not make a difference). I also get a <code>RuntimeWarning: overflow encountered in double_scalars return 1+((10*(A^n))/((A^n)+1))</code> along the way (this could possibly be a division by zero or an infinity) but I am not sure where it is crashing and what from. I've tried changing the parameter ranges to be very small for the sensitivity analysis but alas nothing removes this error.</p>
  <p>Here is my code:</p>
  <pre><code>from scipy import integrate as sp
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.integrate import solve_ivp
  %matplotlib inline
  fig = plt.figure(figsize=(20,10))
  
  def H(A):
      if A&lt;0:
          return 0
      else:
          return 1+((10*(A**n))/((A**n)+1))
      
  def Hill(A):
      if A&lt;0:
          return 0
      else:
          return (A**n)/((A**n)+1)
  
  def NDMsignals(time,Z,iota,upsilon,vartheta,varphi1,varphi2):
      dZdt = [(H(Z[0])*iota)-(upsilon*Z[0])+(Z[1]*varphi1), #A1
      (H(Z[1])*vartheta)-(upsilon*Z[1])+(Z[0]*varphi2)] #A2
      return dZdt
  
  def main(iota,upsilon,vartheta,varphi1,varphi2):
      alpha1 = 0.1
      alpha2 = 0.1
      gamma = 0.75
      n=2.2
      d1 = 1.0
      d2 = 1.0
      beta1 = 10*alpha1
      beta2 = 10*alpha2
      mu1 = 6.0
      mu2 = 6.0
      Upsilon1 = 1.0
      Upsilon2 = 1.0
      c1not = 1.0
      c2not = 1.0
      b1 = 1.0
      b2 = 1.0
      k1 = 0.4
      k2 = 0.4
      sigma1 = 1.0
      sigma2 = 1.0
      tau1 = 0.75
      tau2 = 0.75
  
      time=np.linspace(0,1000,100)
  
      Zinitial = [0.1, 0.1]
  
  # Dimensional Parameters
  
      u = (gamma+b1+b2)
      psi = (tau1/tau2)
      c1eq = (sigma1*k1)/(mu1-k1)
      c2eq = (sigma2*k2)/(mu2-k2)
      xeq = (Upsilon1/k1)*((d1*(c1not-c1eq))-(d2*c1eq)+(d2*c2eq))
      yeq = (Upsilon2/k2)*((d1*(c2not-c2eq))-(d2*c2eq)+(d2*c1eq))
  
      # NDM Parameters
     
      #iota = (xeq)/(tau1)
      #upsilon = (u/alpha1)
      #vartheta = (alpha2*yeq)/(alpha1*tau2)
      #varphi1 = (b1/(psi*alpha1))
      #varphi2 = ((b1*psi)/alpha1)
      
      sol = solve_ivp(lambda t, Z: NDMsignals(time,Z,iota,upsilon,vartheta,varphi1,varphi2), 
                      [time[0],time[-1]], Zinitial, method='BDF', t_eval=time)
      #print([sol.y[0][-1], sol.y[1][-1]])
      return [sol.y[0][-1], sol.y[1][-1], Hill(sol.y[0][-1]), Hill(sol.y[1][-1])]
  
  # Sensitivity Analysis
  
  from SALib.sample import saltelli
  from SALib.analyze import sobol
  
  problem = {
      'num_vars': 5,
      'names':['iota', 'upsilon', 'vartheta', 'varphi1', 'varphi2'],
      'bounds': [[0.01, 10],
                [0.01, 35],
                [0.01, 10],
                [0.01, 3],
                [0.01, 3]]
  }
  
  param_values = saltelli.sample(problem, 200)
  
  Y1 = np.zeros([param_values.shape[0]])
  Y2 = np.zeros([param_values.shape[0]])
  Y3 = np.zeros([param_values.shape[0]])
  Y4 = np.zeros([param_values.shape[0]])
  
  for i, X in enumerate(param_values):
      xx = main(X[0], X[1], X[2], X[3], X[4])
      Y1[i] = xx[0]
      Y2[i] = xx[1]
      Y3[i] = xx[2]
      Y4[i] = xx[3]
  
  #Si1 = sobol.analyze(problem, Y1, print_to_console=True)
  #Si2 = sobol.analyze(problem, Y2, print_to_console=True)
  #Si3 = sobol.analyze(problem, Y3, print_to_console=True)
  #Si4 = sobol.analyze(problem, Y4, print_to_console=True)
  
  plt.scatter(Y3,Y4)
  </code></pre>
  <p>Any help would be fantastic. Thank you.</p>