code: def numdiff(fun, args):
      
  
      # vectorized version
  
      epsilon = 1e-8
      args = list(args)
      v0 = fun(*args)
      N = v0.shape[0]
      l_v = len(v0)
      dvs = []
      for i, a in enumerate(args):
          l_a = (a).shape[1]
          dv = numpy.zeros((N, l_v, l_a))
          nargs = list(args)  #.copy()
          for j in range(l_a):
              xx = args[i].copy()
              xx[:, j] += epsilon
              nargs[i] = xx
              dv[:, :, j] = (fun(*nargs) - v0) / epsilon
          dvs.append(dv)
      return [v0] + dvs
docs: Vectorized numerical differentiation