diff --git a/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt b/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
index f7b284e89..a01ea7fe2 100644
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
@@ -22,8 +22,8 @@ fun main() {
         return DoubleBuffer(x.size) { i ->
             val h = sigma[i] / 5
             val dVector = DoubleBuffer(x.size) { if (it == i) h else 0.0 }
-            val f1 = invoke(x + dVector / 2)
-            val f0 = invoke(x - dVector / 2)
+            val f1 = this(x + dVector / 2)
+            val f0 = this(x - dVector / 2)
             (f1 - f0) / h
         }
     }
diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
index 535c6b39e..92bf86128 100644
--- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
+++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
@@ -28,7 +28,7 @@ public class CMIntegrator(
         val integrator = integratorBuilder(integrand)
         val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
         val remainingCalls = maxCalls - integrand.calls
-        val range = integrand.getFeature<IntegrationRange<Double>>()?.range
+        val range = integrand.getFeature<IntegrationRange>()?.range
             ?: error("Integration range is not provided")
         val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive)
 
diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt
index 071bac315..1c9915563 100644
--- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt
+++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt
@@ -17,7 +17,7 @@ public class GaussRuleIntegrator(
 ) : UnivariateIntegrator<Double> {
 
     override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
-        val range = integrand.getFeature<IntegrationRange<Double>>()?.range
+        val range = integrand.getFeature<IntegrationRange>()?.range
             ?: error("Integration range is not provided")
         val integrator: GaussIntegrator = getIntegrator(range)
         //TODO check performance
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt
index bc23e2f1b..749f3cc72 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt
@@ -8,6 +8,13 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Field
 import space.kscience.kmath.structures.*
 
+/**
+ * Set of univariate integration ranges. First components correspond to ranges themselves, second components to number of
+ * integration nodes per range
+ */
+public class UnivariateIntegrandRanges(public val ranges: List<Pair<ClosedRange<Double>, Int>>) : IntegrandFeature {
+    public constructor(vararg pairs: Pair<ClosedRange<Double>, Int>) : this(pairs.toList())
+}
 
 /**
  * A simple one-pass integrator based on Gauss rule
@@ -18,9 +25,29 @@ public class GaussIntegrator<T : Any>(
 
     private fun buildRule(integrand: UnivariateIntegrand<T>): Pair<Buffer<Double>, Buffer<Double>> {
         val factory = integrand.getFeature<GaussIntegratorRuleFactory>() ?: GaussLegendreRuleFactory
-        val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
-        val range = integrand.getFeature<IntegrationRange<Double>>()?.range ?: 0.0..1.0
-        return factory.build(numPoints, range)
+        val predefinedRanges = integrand.getFeature<UnivariateIntegrandRanges>()
+        if (predefinedRanges == null || predefinedRanges.ranges.isEmpty()) {
+            val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+            val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+            return factory.build(numPoints, range)
+        } else {
+            val ranges = predefinedRanges.ranges
+            var counter = 0
+            val length = ranges.sumOf { it.second }
+            val pointsArray = DoubleArray(length)
+            val weightsArray = DoubleArray(length)
+
+            for (range in ranges) {
+                val rule = factory.build(range.second, range.first)
+                repeat(rule.first.size) { i ->
+                    pointsArray[counter] = rule.first[i]
+                    weightsArray[counter] = rule.second[i]
+                    counter++
+                }
+
+            }
+            return pointsArray.asBuffer() to weightsArray.asBuffer()
+        }
     }
 
     override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
@@ -49,7 +76,8 @@ public class GaussIntegrator<T : Any>(
  * Following features are evaluated:
  * * [GaussIntegratorRuleFactory] - A factory for computing the Gauss integration rule. By default uses [GaussLegendreRuleFactory]
  * * [IntegrationRange] - the univariate range of integration. By default uses 0..1 interval.
- * * [IntegrandMaxCalls] - the maximum number of function calls during integration. For non-iterative rules, always uses the maximum number of points. By default uses 100 points.
+ * * [IntegrandMaxCalls] - the maximum number of function calls during integration. For non-iterative rules, always uses the maximum number of points. By default uses 10 points.
+ * * [UnivariateIntegrandRanges] - Set of ranges and number of points per range. Defaults to given [IntegrationRange] and [IntegrandMaxCalls]
  */
 @UnstableKMathAPI
 public fun <T : Any> Field<T>.integrate(
@@ -64,15 +92,25 @@ public fun <T : Any> Field<T>.integrate(
 @UnstableKMathAPI
 public fun <T : Any> Field<T>.integrate(
     range: ClosedRange<Double>,
-    numPoints: Int = 100,
+    order: Int = 10,
+    intervals: Int = 10,
     vararg features: IntegrandFeature,
     function: (Double) -> T,
-): UnivariateIntegrand<T> = GaussIntegrator(this).integrate(
-    UnivariateIntegrand(
-        function,
-        IntegrationRange(range),
-        GaussLegendreRuleFactory,
-        IntegrandMaxCalls(numPoints),
-        *features
+): UnivariateIntegrand<T> {
+    require(range.endInclusive > range.start) { "The range upper bound should be higher than lower bound" }
+    require(order > 1) { "The order of polynomial must be more than 1" }
+    require(intervals > 0) { "Number of intervals must be positive" }
+    val rangeSize = (range.endInclusive - range.start) / intervals
+    val ranges = UnivariateIntegrandRanges(
+        (0 until intervals).map { i -> (rangeSize * i)..(rangeSize * i + 1) to order }
     )
-)
\ No newline at end of file
+    return GaussIntegrator(this).integrate(
+        UnivariateIntegrand(
+            function,
+            IntegrationRange(range),
+            GaussLegendreRuleFactory,
+            ranges,
+            *features
+        )
+    )
+}
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
index 0b41a3f8b..bcd5005c4 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
@@ -33,7 +33,7 @@ public fun <T : Any> UnivariateIntegrand(
 public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
 
 @JvmInline
-public value class IntegrationRange<T : Comparable<T>>(public val range: ClosedRange<T>) : IntegrandFeature
+public value class IntegrationRange(public val range: ClosedRange<Double>) : IntegrandFeature
 
 public val <T : Any> UnivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value