diff --git a/CHANGELOG.md b/CHANGELOG.md
index 77ee2909e..7c093ade1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,12 +13,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 - Add NURBS surface to PPM conversion examples.
 - Utilized ForUnitTest for testing.
 - Added Support for `OpenMP` and `do concurrent`.
+- Implemented memory cleanup in the examples and tests.
 
 ### Changed
 
 - Updated `README.md` file.
 - Updated tests to use ForUnitTest.
 - Added ParaView to the list of References in the README file.
+- Used `matmul` instead of `dot_product` in the `put_to_nurbs` subroutine.
 
 ### Removed
 
diff --git a/README.md b/README.md
index b9d57743b..613af01af 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@
 
 ![](logo/logo.png)
 
-**ForCAD**: A Parallel Fortran library for Geometric Modeling using NURBS (Non-Uniform Rational B-Splines).
+**ForCAD**: A parallel Fortran library for geometric modeling using NURBS (Non-Uniform Rational B-Splines).
 
 ForCAD supports **B-Spline**, **NURBS**, **Bezier**, and **Rational Bezier** curves, surfaces, and volumes.
 
diff --git a/example/demo_curve.f90 b/example/demo_curve.f90
index 45cf2f8c1..b6b4ec08e 100644
--- a/example/demo_curve.f90
+++ b/example/demo_curve.f90
@@ -23,6 +23,9 @@ program example_nurbs_curve
     !> Set control points and weights for the NURBS curve object
     call nurbs%set(Xc, Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export initial control points to a VTK file
     call nurbs%export_Xc('vtk/demo_curve_Xc.vtk')
 
diff --git a/example/demo_surface.f90 b/example/demo_surface.f90
index 7f1a8403e..118ceacd1 100644
--- a/example/demo_surface.f90
+++ b/example/demo_surface.f90
@@ -23,6 +23,9 @@ program example_nurbs_surface
     !> Set control points and weights for the NURBS surface object
     call nurbs%set([10,10],Xc,Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export initial control points to a VTK file
     call nurbs%export_Xc('vtk/demo_surface_Xc.vtk')
 
diff --git a/example/demo_volume.f90 b/example/demo_volume.f90
index 84dabe27b..86b9e4511 100644
--- a/example/demo_volume.f90
+++ b/example/demo_volume.f90
@@ -23,6 +23,9 @@ program example_nurbs_volume
     !> Set control points and weights for the NURBS volume object
     call nurbs%set([2,2,2], Xc, Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export initial control points to a VTK file
     call nurbs%export_Xc('vtk/demo_volume_Xc.vtk')
 
diff --git a/example/example_curve_1.f90 b/example/example_curve_1.f90
index a0bcfe2eb..9a57986a6 100644
--- a/example/example_curve_1.f90
+++ b/example/example_curve_1.f90
@@ -31,6 +31,9 @@ program example1_curve
     !> Wc is optional
     call nurbs%set(knot, Xc, Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export control points to a VTK file
     call nurbs%export_Xc('vtk/nurbs_curve_Xc.vtk')
 
diff --git a/example/example_ppm1.f90 b/example/example_ppm1.f90
index 736e10262..8817bc6db 100644
--- a/example/example_ppm1.f90
+++ b/example/example_ppm1.f90
@@ -354,4 +354,14 @@ program example_ppm1
     call image%finalize()
     call t%timer_stop(message='Saving the image')
 
+
+
+
+
+
+
+    ! Clean up
+    ! call cmap%finalize()
+    deallocate(px, Xg, z_values)
+
 end program
diff --git a/example/example_ppm2.f90 b/example/example_ppm2.f90
index c4c7f700a..4ccfc38be 100644
--- a/example/example_ppm2.f90
+++ b/example/example_ppm2.f90
@@ -363,4 +363,14 @@ program example_ppm2
     call image%finalize()
     call t%timer_stop(message='Saving the image')
 
+
+
+
+
+
+
+    ! Clean up
+    ! call cmap%finalize()
+    deallocate(px, Xg, z_values)
+
 end program
diff --git a/example/example_ppm3.f90 b/example/example_ppm3.f90
index 289ed9458..51696f45e 100644
--- a/example/example_ppm3.f90
+++ b/example/example_ppm3.f90
@@ -346,4 +346,14 @@ program example_ppm3
     call image%finalize()
     call t%timer_stop(message='Saving the image')
 
+
+
+
+
+
+
+    ! Clean up
+    ! call cmap%finalize()
+    deallocate(px, Xg, z_values)
+
 end program
diff --git a/example/example_surface_1.f90 b/example/example_surface_1.f90
index 10b754dab..31db38343 100644
--- a/example/example_surface_1.f90
+++ b/example/example_surface_1.f90
@@ -30,6 +30,9 @@ program example3_surface
     !> Set knot vectors, control points, and weights for the NURBS surface object
     call nurbs%set(knot1, knot2, Xc, Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export the control points to a VTK file
     call nurbs%export_Xc('vtk/nurbs_surface_Xc.vtk')
 
diff --git a/example/example_volume_1.f90 b/example/example_volume_1.f90
index d31f6f99b..39fef0afc 100644
--- a/example/example_volume_1.f90
+++ b/example/example_volume_1.f90
@@ -30,6 +30,9 @@ program example3_volume
     !> Wc is optional.
     call nurbs%set(knot1, knot2, knot3, Xc, Wc)
 
+    !> Deallocate local arrays
+    deallocate(Xc, Wc)
+
     !> Export the control points to a VTK file
     call nurbs%export_Xc('vtk/nurbs_volume_Xc.vtk')
 
diff --git a/example/put_to_nurbs.f90 b/example/put_to_nurbs.f90
index 829bce711..c68b7e590 100644
--- a/example/put_to_nurbs.f90
+++ b/example/put_to_nurbs.f90
@@ -35,6 +35,9 @@ program example_morph
     !> Map the shape into the shape
     call control_shape%put_to_nurbs(X, elem)
 
+    !> Deallocate local variables
+    deallocate(X, elem)
+
     !> Export the shape and the control shape to vtk files
     call control_shape%export_Xc('vtk/control_shape.vtk')
     call control_shape%export_Xg('vtk/shape.vtk')
@@ -47,4 +50,7 @@ program example_morph
     !> Show the control geometry and geometry using PyVista
     call control_shape%show('vtk/control_shape.vtk','vtk/shape.vtk')
 
+    !> Finalize the control shape
+    call control_shape%finalize()
+
 end program
\ No newline at end of file
diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90
index fd983bf69..a08ce2187 100644
--- a/src/forcad_nurbs_surface.f90
+++ b/src/forcad_nurbs_surface.f90
@@ -2103,6 +2103,8 @@ impure function compute_Xg_nurbs_2d(Xt, knot1, knot2, degree, nc, ng, Xc, Wc) re
     integer :: i
 
     allocate(Xg(ng(1)*ng(2), size(Xc,2)))
+    allocate(Tgc(nc(1)*nc(2)))
+
     !$OMP PARALLEL DO PRIVATE(Tgc)
     do i = 1, ng(1)*ng(2)
         Tgc = kron(&
diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90
index 82380a42e..509870ba2 100644
--- a/src/forcad_nurbs_volume.f90
+++ b/src/forcad_nurbs_volume.f90
@@ -294,7 +294,8 @@ pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f
         end if
 
         if (allocated(this%Xg)) deallocate(this%Xg)
-
+        allocate(this%Xg(this%ng(1)*this%ng(2)*this%ng(3), size(this%Xc,2)))
+        
         if (allocated(this%Wc)) then ! NURBS volume
             this%Xg = compute_Xg_nurbs_3d(Xt_, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc, this%Wc)
         else
@@ -1813,6 +1814,7 @@ pure subroutine put_to_nurbs(this, X, elemConn)
         Xt(:,3) = (X(:,3) -  min_X3) / (max_X3 -  min_X3)
 
         allocate(this%Xg(size(Xt,1), size(this%Xc,2)))
+        allocate(Tgc(this%nc(1)*this%nc(2)*this%nc(3)))
 
         if (allocated(this%Wc)) then ! NURBS volume
             do i = 1, size(Xt, 1)
@@ -1821,9 +1823,7 @@ pure subroutine put_to_nurbs(this, X, elemConn)
                 Tgc3 = basis_bspline(Xt(i,3), this%knot3, this%nc(3), this%degree(3))
                 Tgc = kron(Tgc3, kron(Tgc2, Tgc1))
                 Tgc = Tgc*(this%Wc/(dot_product(Tgc,this%Wc)))
-                do j = 1, size(this%Xc, 2)
-                    this%Xg(i,j) = dot_product(Tgc,this%Xc(:,j))
-                end do
+                this%Xg(i,:) = matmul(Tgc,this%Xc)
             end do
         else ! B-Spline volume
             do i = 1, size(Xt, 1)
@@ -1831,9 +1831,7 @@ pure subroutine put_to_nurbs(this, X, elemConn)
                 Tgc2 = basis_bspline(Xt(i,2), this%knot2, this%nc(2), this%degree(2))
                 Tgc3 = basis_bspline(Xt(i,3), this%knot3, this%nc(3), this%degree(3))
                 Tgc = kron(Tgc3, kron(Tgc2, Tgc1))
-                do j = 1, size(this%Xc, 2)
-                    this%Xg(i,j) = dot_product(Tgc,this%Xc(:,j))
-                end do
+                this%Xg(i,:) = matmul(Tgc,this%Xc)
             end do
         end if
 
@@ -2596,6 +2594,8 @@ impure function compute_Xg_nurbs_3d(Xt, knot1, knot2, knot3, degree, nc, ng, Xc,
     integer :: i
 
     allocate(Xg(ng(1)*ng(2)*ng(3), size(Xc,2)))
+    allocate(Tgc(nc(1)*nc(2)*nc(3)))
+
     !$OMP PARALLEL DO PRIVATE(Tgc)
     do i = 1, ng(1)*ng(2)*ng(3)
         Tgc = kron(basis_bspline(Xt(i,3), knot3, nc(3), degree(3)),&
diff --git a/src/forcad_utils.f90 b/src/forcad_utils.f90
index a66734501..7e742166a 100644
--- a/src/forcad_utils.f90
+++ b/src/forcad_utils.f90
@@ -248,7 +248,7 @@ pure function repelem(a, b) result(c)
             c(l+1:l+n) = a(i)
             l = l + n
         end do
-    end function repelem
+    end function
     !===============================================================================
 
 
@@ -328,18 +328,16 @@ pure subroutine cmp_elemConn_Cn_L(nnode, p, Xth,vecKnot_mul, elemConn)
         real(rk), intent(in), contiguous :: Xth(:)
         integer, allocatable, intent(out) :: elemConn(:,:)
         integer, allocatable :: nodes(:)
-        integer :: i, l, nnel, m, nelem
+        integer :: i, nnel, m, nelem
 
         nnel = p + 1
         nodes = [(i, i=1, nnode)]
         nelem = size(Xth) - 1
         allocate(elemConn(nelem,nnel))
-        l = 0
         m = -p
         do i = 1, nelem
             m = m + vecKnot_mul(i)
-            l = l + 1
-            elemConn(l,:) = nodes(m:m+p)
+            elemConn(i,:) = nodes(m:m+p)
         end do
     end subroutine
     !===============================================================================
diff --git a/test/test_curve.f90 b/test/test_curve.f90
index 193e37451..97a7340f7 100644
--- a/test/test_curve.f90
+++ b/test/test_curve.f90
@@ -58,5 +58,6 @@ program test_curve
 
     call nurbs%finalize()
     call bsp%finalize()
+    deallocate(Xc, Wc, Xg, Xgb)
 
 end program
diff --git a/test/test_surface.f90 b/test/test_surface.f90
index abe8718a7..c593404cb 100644
--- a/test/test_surface.f90
+++ b/test/test_surface.f90
@@ -66,6 +66,7 @@ program test_surface
 
     call nurbs%finalize()
     call bsp%finalize()
+    deallocate(Xc, Wc, Xg, Xgb)
 
 contains
 
diff --git a/test/test_volume.f90 b/test/test_volume.f90
index 81aa2586d..c7eb8e83d 100644
--- a/test/test_volume.f90
+++ b/test/test_volume.f90
@@ -72,6 +72,7 @@ program test_volume
 
     call nurbs%finalize()
     call bsp%finalize()
+    deallocate(Xc, Wc, Xg, Xgb)
 
 contains