diff --git a/example/demo.f90 b/example/demo.f90
index 1ada4e2..89a09d6 100644
--- a/example/demo.f90
+++ b/example/demo.f90
@@ -21,15 +21,16 @@
 ! SOFTWARE.
 !-------------------------------------------------------------------------------
 ! Contributed by vmagnin: 2023-09-26
-! Last modification: vmagnin 2023-10-19
+! Last modification: gha3mi 2024-01-28
 !-------------------------------------------------------------------------------
 
 program demo
-    use forcolormap, only: Colormap, colormaps_list, wp
+    use forcolormap, only: Colormap, colormaps_list, wp, bezier
     implicit none
 
     integer :: i
     type(Colormap) :: cmap, custom_cmap
+    integer, allocatable :: colors(:,:)
     ! A discrete colormap with 8 levels, from black to white:
     integer, dimension(0:7, 3) :: my_colormap = reshape( [ &
           0,     0,     0,   &
@@ -68,6 +69,16 @@ program demo
     call test_colormap(custom_cmap, 'a_loaded_colormap_test', encoding='binary')
     call custom_cmap%print()
 
+    ! You can also create your own colormap using bezier interpolation:
+    ! Define control colors.
+    allocate(colors(3,3))
+    colors(1,:) = [255, 0, 0] ! Red
+    colors(2,:) = [0, 255, 0] ! Green
+    colors(3,:) = [0, 0, 255] ! Blue
+    call custom_cmap%create('custom_bezier', 0.0_wp, 2.0_wp, bezier(colors, levels=1024)) ! levels is optional, default is 256
+    call custom_cmap%colorbar('custom_colorbar_bezier', encoding='binary')
+    call test_colormap(custom_cmap, 'custom_test_bezier', encoding='binary') 
+
     contains
 
     subroutine test_colormap(self, filename, encoding)
diff --git a/src/colormap_class.f90 b/src/colormap_class.f90
index ea1c764..040445e 100644
--- a/src/colormap_class.f90
+++ b/src/colormap_class.f90
@@ -21,7 +21,7 @@
 ! SOFTWARE.
 !-------------------------------------------------------------------------------
 ! Contributed by vmagnin: 2023-09-26
-! Last modification: gha3mi 2024-01-06
+! Last modification: gha3mi 2024-01-28
 !-------------------------------------------------------------------------------
 
 
@@ -34,7 +34,7 @@ module forcolormap
     implicit none
     private
 
-    public :: wp
+    public :: wp, bezier
 
     ! List of built-in colormaps:
     character(*), dimension(*), public, parameter :: colormaps_list = &
@@ -927,4 +927,44 @@ pure subroutine extract(self, extractedLevels, name, zmin, zmax, reverse)
         end if
     end subroutine extract
 
+    ! Create colormap from continuous Bezier interpolation of control colors
+    pure function bezier(colors, levels) result(map)
+        integer, dimension(:,:), intent(in) :: colors
+        integer, intent(in), optional :: levels
+        integer, dimension(:,:), allocatable :: map
+        real(wp), dimension(:,:), allocatable :: map_r
+        integer :: order, i, j, levels_
+        real(wp) :: t
+
+        ! Set default value for levels
+        if (present(levels)) then
+            levels_ = levels
+        else
+            levels_ = 256
+        end if
+
+        ! Order of the Bezier curve
+        order = size(colors, 1) - 1
+
+        allocate(map_r(levels_,3), map(levels_,3)) ! 3 for RGB
+        do i = 1,levels_
+            t = real(i-1, wp) / real(levels_-1, wp)
+            map_r(i,:) = 0.0_wp
+            do j = 0, order
+                map_r(i,:) = map_r(i,:) + real(colors(j+1,:), wp)*&
+                real(factorial(order), wp)/(real(factorial(j), wp)*real(factorial(order-j), wp)) * t**j * (1.0_wp-t)**(order-j)
+            end do
+            map(i,:) = scale_real_int(map_r(i,:), 0, 255) ! Scale to integer RGB range
+        end do
+    end function bezier
+
+    ! Factorial function used for Bezier interpolation
+    pure function factorial(n) result(result)
+        integer, intent(in) :: n
+        integer :: result, i
+        result = 1
+        do concurrent (i = 2:n)
+            result = result * i
+        end do
+    end function factorial
 end module forcolormap