diff --git a/example/linalg/example_repmat.f90 b/example/linalg/example_repmat.f90
new file mode 100644
index 000000000..c907ec879
--- /dev/null
+++ b/example/linalg/example_repmat.f90
@@ -0,0 +1,12 @@
+program example_repmat
+  use stdlib_linalg, only: repmat
+  implicit none
+  real, allocatable :: a(:, :),b(:,:)
+  a = reshape([1., 2., 3., 4.],[2, 2],order=[2, 1])
+  b = repmat(a, 2, 2)
+! A = reshape([1., 2., 1., 2.,&
+!              3., 4., 3., 4.,&
+!              1., 2., 1., 2.,&
+!              3., 4., 3., 4.,&
+!               ], [4, 4],order=[2, 1])
+end program example_repmat
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0fb95a2d3..a056e791b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -22,6 +22,7 @@ set(fppFiles
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp
     stdlib_linalg_outer_product.fypp
+    stdlib_linalg_repmat.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index bc1017f0a..d064c948e 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -14,6 +14,7 @@ module stdlib_linalg
   public :: eye
   public :: trace
   public :: outer_product
+  public :: repmat
   public :: is_square
   public :: is_diagonal
   public :: is_symmetric
@@ -92,6 +93,21 @@ module stdlib_linalg
     #:endfor
   end interface outer_product
 
+  ! repeat matrix (of 2d array)
+  interface repmat
+    !! version: experimental
+    !!
+    !! Creates large matrices from a small array, `repmat()` repeats the given values of the array to create the large matrix.
+    !! ([Specification](../page/specs/stdlib_linalg.html#
+    !! repmat-creates-large-matrices-from-a-small-array))
+    #:for k1, t1 in RCI_KINDS_TYPES
+      pure module function repmat_${t1[0]}$${k1}$(a,m,n) result(res)
+        ${t1}$, intent(in) :: a(:,:)
+        integer, intent(in) :: m,n
+        ${t1}$ :: res(m*size(a,1),n*size(a,2))
+      end function repmat_${t1[0]}$${k1}$
+    #:endfor
+  end interface repmat
 
   ! Check for squareness
   interface is_square
diff --git a/src/stdlib_linalg_repmat.fypp b/src/stdlib_linalg_repmat.fypp
new file mode 100644
index 000000000..3ea47a064
--- /dev/null
+++ b/src/stdlib_linalg_repmat.fypp
@@ -0,0 +1,30 @@
+#:include "common.fypp"
+#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_repmat
+
+  implicit none
+
+contains
+
+  #:for k1, t1 in RCI_KINDS_TYPES
+    pure module function repmat_${t1[0]}$${k1}$(a, m, n) result(res)
+      ${t1}$, intent(in) :: a(:,:)
+      integer, intent(in) :: m,n
+      ${t1}$ :: res(m*size(a,1),n*size(a,2))
+      integer ::i,j,k,l
+      associate(ma=>size(a,1),na=>size(a,2))
+         do j=1,n
+            do l=1,na
+               do i=1,m
+                  do k=1,ma
+                     res((i-1)*ma+k,(j-1)*na+l)=a(k,l)
+                  end do
+               end do
+            end do
+         end do
+      end associate
+    end function repmat_${t1[0]}$${k1}$
+
+  #:endfor
+
+end submodule
diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp
index f74cbff6b..f278345da 100644
--- a/test/linalg/test_linalg.fypp
+++ b/test/linalg/test_linalg.fypp
@@ -3,7 +3,7 @@
 module test_linalg
     use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
     use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
-    use stdlib_linalg, only: diag, eye, trace, outer_product
+    use stdlib_linalg, only: diag, eye, trace, outer_product, repmat
 
     implicit none
 
@@ -57,7 +57,17 @@ contains
             new_unittest("outer_product_int8", test_outer_product_int8), &
             new_unittest("outer_product_int16", test_outer_product_int16), &
             new_unittest("outer_product_int32", test_outer_product_int32), &
-            new_unittest("outer_product_int64", test_outer_product_int64) &
+            new_unittest("outer_product_int64", test_outer_product_int64), &
+            new_unittest("test_repmat_rsp",test_repmat_rsp), &
+            new_unittest("test_repmat_rdp",test_repmat_rdp), &
+            new_unittest("test_repmat_rqp",test_repmat_rqp), &
+            new_unittest("test_repmat_csp",test_repmat_csp), &
+            new_unittest("test_repmat_cdp",test_repmat_cdp), &
+            new_unittest("test_repmat_cqp",test_repmat_cqp), &
+            new_unittest("test_repmat_int8",test_repmat_int8), &
+            new_unittest("test_repmat_int16",test_repmat_int16), &
+            new_unittest("test_repmat_int32",test_repmat_int32), &
+            new_unittest("test_repmat_int64",test_repmat_int64) &
             ]
 
     end subroutine collect_linalg
@@ -703,6 +713,197 @@ contains
     end subroutine test_outer_product_int64
 
 
+    subroutine test_repmat_rqp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+#:if WITH_QP
+       integer, parameter :: n = 2
+       real(qp) :: a(n,n)
+       real(qp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+#:else
+        call skip_test(error, "Quadruple precision is not enabled")
+#:endif
+    end subroutine test_repmat_rqp
+    subroutine test_repmat_cqp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+#:if WITH_QP
+       integer, parameter :: n = 2
+       complex(qp) :: a(n,n)
+       complex(qp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+#:else
+        call skip_test(error, "Quadruple precision is not enabled")
+#:endif
+    end subroutine test_repmat_cqp
+
+
+    subroutine test_repmat_rsp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       real(sp) :: a(n,n)
+       real(sp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_rsp
+
+    subroutine test_repmat_rdp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       real(dp) :: a(n,n)
+       real(dp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_rdp
+
+
+    subroutine test_repmat_csp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       complex(sp) :: a(n,n)
+       complex(sp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_csp
+
+    subroutine test_repmat_cdp(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       complex(dp) :: a(n,n)
+       complex(dp),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_cdp
+   
+    subroutine test_repmat_int8(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       integer(int8) :: a(n,n)
+       integer(int8),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_int8
+
+    subroutine test_repmat_int16(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       integer(int16) :: a(n,n)
+       integer(int16),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_int16
+
+    subroutine test_repmat_int32(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       integer(int32) :: a(n,n)
+       integer(int32),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_int32
+
+    subroutine test_repmat_int64(error)
+       !> Error handling
+       type(error_type), allocatable, intent(out) :: error
+       integer, parameter :: n = 2
+       integer(int64) :: a(n,n)
+       integer(int64),allocatable:: expected(:,:),diff(:,:)
+
+       a=reshape([1,2,3,4],shape(a),order=[2,1])
+       expected = reshape([&
+          1,2,1,2,&
+          3,4,3,4,&
+          1,2,1,2,&
+          3,4,3,4 ],[4,4],order=[2,1])
+       diff = expected - repmat(a, 2, 2)
+       call check(error, all(abs(diff) == 0), &
+          "all(abs(diff) == 0) failed.")
+    end subroutine test_repmat_int64
+
+
+
     pure recursive function catalan_number(n) result(value)
         integer, intent(in) :: n
         integer :: value