diff --git a/integration_tests/CMakeLists.txt b/integration_tests/CMakeLists.txt
index a6ed37cd1b..f38ea7fb0b 100644
--- a/integration_tests/CMakeLists.txt
+++ b/integration_tests/CMakeLists.txt
@@ -409,6 +409,8 @@ RUN(NAME test_str_comparison LABELS cpython llvm c)
 RUN(NAME test_bit_length     LABELS cpython llvm c)
 RUN(NAME str_to_list_cast    LABELS cpython llvm c)
 
+RUN(NAME test_package_01     LABELS cpython llvm)
+
 RUN(NAME generics_01         LABELS cpython llvm c)
 RUN(NAME generics_02         LABELS cpython llvm c)
 RUN(NAME generics_array_01   LABELS cpython llvm c)
diff --git a/integration_tests/nrp/__init__.py b/integration_tests/nrp/__init__.py
new file mode 100644
index 0000000000..f5631d2af1
--- /dev/null
+++ b/integration_tests/nrp/__init__.py
@@ -0,0 +1 @@
+from .nr import newton_raphson
diff --git a/integration_tests/nrp/nr.py b/integration_tests/nrp/nr.py
new file mode 100644
index 0000000000..89be2fb0d0
--- /dev/null
+++ b/integration_tests/nrp/nr.py
@@ -0,0 +1,20 @@
+from ltypes import i32, f64
+
+
+def func(x: f64, c: f64) -> f64:
+    return x**2.0 - c**2.0
+
+
+def func_prime(x: f64) -> f64:
+    return 2.0*x
+
+
+def newton_raphson(x: f64, c: f64, maxiter: i32) -> f64:
+    h: f64
+    err: f64 = 1e-5
+    i: i32 = 0
+    while abs(func(x, c)) > err and i < maxiter:
+        h = func(x, c) /  func_prime(x)
+        x -= h
+        i += 1
+    return x
diff --git a/integration_tests/test_package_01.py b/integration_tests/test_package_01.py
new file mode 100644
index 0000000000..0026a0c8c7
--- /dev/null
+++ b/integration_tests/test_package_01.py
@@ -0,0 +1,13 @@
+from nrp import newton_raphson
+from ltypes import f64, i32
+
+
+def check():
+    x0: f64 = 20.0
+    c: f64 = 3.0
+    maxiter: i32 = 20
+    x: f64
+    x = newton_raphson(x0, c, maxiter)
+    assert abs(x - 3.0) < 1e-5
+
+check()