code: def hinit(func, x, t, pos_neg, f0, iord, hmax, rtol, atol, args):
      
      sk = atol + rtol * np.fabs(x)
      dnf = np.sum(np.square(f0 / sk), axis=0)
      dny = np.sum(np.square(x / sk), axis=0)
  
      h = np.sqrt(dny / dnf) * 0.01
  
      h = np.min([h, np.fabs(hmax)])
      h = custom_sign(h, pos_neg)
  
      # perform an explicit Euler step
      xx1 = x + h * f0
      f1 = np.array(func(xx1, t[0] + h, *args))
  
      # estimate the second derivative of the solution
      der2 = np.sum(np.square((f1 - f0) / sk), axis=0)
      der2 = np.sqrt(der2) / h
  
      # step size is computed such that h ** iord * max_d(norm(f0), norm(der2)) = 0.01
      der12 = np.max([np.fabs(der2), np.sqrt(dnf)])
      h1 = np.power(0.01 / der12, 1.0 / iord)
  
      h = np.min([100.0 * np.fabs(h), np.min([np.fabs(h1), np.fabs(hmax)])])
  
      return custom_sign(h, pos_neg), f0, f1, xx1
docs: Estimate initial step size