diff --git a/source/module_base/math_polyint.h b/source/module_base/math_polyint.h
index 0cc13f8f83..b15f50678c 100644
--- a/source/module_base/math_polyint.h
+++ b/source/module_base/math_polyint.h
@@ -18,6 +18,19 @@ class PolyInt
     //========================================================
     // Polynomial_Interpolation
     //========================================================
+
+    /**
+     * @brief Lagrange interpolation
+     * 
+     * @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
+     * @param dim1 [in] index of 1st dimension of table/y
+     * @param dim2 [in] index of 2nd dimension of table/y
+     * @param y [out] three dimension matrix to store the predicted value
+     * @param dim_y [in] index of 3rd dimension of y to store predicted value
+     * @param table_length [in] length of 3rd dimension of table
+     * @param table_interval [in] interval of 3rd dimension of table
+     * @param x [in] the position in 3rd dimension to be predicted
+     */
     static void Polynomial_Interpolation
     (
         const ModuleBase::realArray &table,
@@ -30,6 +43,17 @@ class PolyInt
         const double &x
     );
 
+    /**
+     * @brief Lagrange interpolation
+     * 
+     * @param table [in] three dimension matrix, the data in 3rd dimension is used to do prediction
+     * @param dim1 [in] index of 1st dimension of table
+     * @param dim2 [in] index of 2nd dimension of table
+     * @param table_length [in] length of 3rd dimension of table
+     * @param table_interval [in] interval of 3rd dimension of table
+     * @param x [in] the position in 3rd dimension to be predicted
+     * @return double the predicted value
+     */
     static double Polynomial_Interpolation
     (
         const ModuleBase::realArray &table,
@@ -37,10 +61,24 @@ class PolyInt
         const int &dim2,
         const int &table_length,
         const double &table_interval,
-        const double &x             // input value
+        const double &x            
     );
 
-    static double Polynomial_Interpolation             // pengfei Li 2018-3-23
+    /**
+     * @brief Lagrange interpolation
+     * 
+     * @param table [in] four dimension matrix, the data in 4th dimension is used to do prediction
+     * @param dim1 [in] index of 1st dimension of table
+     * @param dim2 [in] index of 2nd dimension of table 
+     * @param dim3 [in] index of 3rd dimension of table 
+     * @param table_length [in] length of 4th dimension of table
+     * @param table_interval [in] interval of 4th dimension of table
+     * @param x [in] the position in 4th dimension to be predicted
+     * @return double the predicted value
+     * @author pengfei Li
+     * @date 2018-3-23
+     */
+    static double Polynomial_Interpolation            
     (
         const ModuleBase::realArray &table,
         const int &dim1,
@@ -48,23 +86,41 @@ class PolyInt
         const int &dim3,
         const int &table_length,
         const double &table_interval,
-        const double &x             // input value
+        const double &x            
     );
 
+    /**
+     * @brief  Lagrange interpolation
+     * 
+     * @param table [in] the data used to do prediction
+     * @param table_length [in] length of table
+     * @param table_interval [in] interval of table
+     * @param x [in] the position to be predicted
+     * @return double the predicted value
+     */
 	static double Polynomial_Interpolation
 	(
         const double *table,
         const int &table_length,
         const double &table_interval,
-        const double &x             // input value
+        const double &x            
     );
 
+    /**
+     * @brief Lagrange interpolation
+     * 
+     * @param xpoint [in] array of postion
+     * @param ypoint [in] array of data to do prediction
+     * @param table_length [in] length of xpoint
+     * @param x [in] position to be predicted
+     * @return double predicted value
+     */
     static double Polynomial_Interpolation_xy
     (
         const double *xpoint,
         const double *ypoint,
         const int table_length,
-        const double &x             // input value
+        const double &x            
     );
 
 };
diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt
index ad36d6940b..c7cbd3c438 100644
--- a/source/module_base/test/CMakeLists.txt
+++ b/source/module_base/test/CMakeLists.txt
@@ -67,4 +67,8 @@ AddTest(
   TARGET base_mathzone
   LIBS ${math_libs}
   SOURCES mathzone_test.cpp ../mathzone.cpp ../matrix3.cpp ../matrix.cpp ../tool_quit.cpp ../global_variable.cpp ../global_file.cpp ../memory.cpp ../timer.cpp
+)
+AddTest(
+  TARGET base_math_polyint
+  SOURCES math_polyint_test.cpp ../math_polyint.cpp ../realarray.cpp ../timer.cpp
 )
\ No newline at end of file
diff --git a/source/module_base/test/math_polyint_test.cpp b/source/module_base/test/math_polyint_test.cpp
new file mode 100644
index 0000000000..bb5a5937e9
--- /dev/null
+++ b/source/module_base/test/math_polyint_test.cpp
@@ -0,0 +1,136 @@
+#include"../math_polyint.h"
+#include"gtest/gtest.h"
+#include"../realarray.h"
+#include<math.h>
+
+#define doublethreshold 1e-9
+
+/************************************************
+*  unit test of class PolyInt
+***********************************************/
+
+/**
+ * This unit test is to verify the accuracy of 
+ * interpolation method on the function sin(x)/x
+ * with a interval of 0.01.
+ * sin(x)/x is one of the solution of spherical bessel
+ * function when l=0.
+ * 
+ * - Tested function:
+ *      - 4 types of Polynomial_Interpolation 
+ *      - Polynomial_Interpolation_xy
+ */
+
+
+class bessell0 : public testing::Test
+{
+    protected:
+    
+    int TableLength = 400;
+    double interval = 0.01;
+    ModuleBase::realArray table3,table4;
+    ModuleBase::realArray y3;
+    double *tablex;
+    double *tabley;
+
+    double sinc(double x) {return sin(x)/x;}
+
+    void SetUp()
+    {
+        tablex = new double[TableLength];
+        tabley = new double[TableLength];
+        table3.create(1,1,TableLength);
+        table4.create(1,1,1,TableLength);
+        y3.create(1,1,TableLength);
+
+        for(int i=1;i<TableLength;++i) 
+        {
+            table3(0,0,i) = sinc(i * interval);
+            table4(0,0,0,i) = sinc(i * interval); 
+            tablex[i] = i * interval;
+            tabley[i] = sinc(i * interval);
+        }
+    }
+
+    void TearDown()
+    {
+        delete [] tablex;
+        delete [] tabley;
+    }
+};
+
+
+TEST_F(bessell0,PolynomialInterpolationThreeDimensionY)
+{
+    ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,1,TableLength,interval,0.1);
+    ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,2,TableLength,interval,1.005);
+    ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,3,TableLength,interval,2.005);
+    ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,4,TableLength,interval,3.005);
+    ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,y3,5,TableLength,interval,3.505);
+
+    EXPECT_NEAR(y3(0,0,1),sinc(0.1),doublethreshold);
+    EXPECT_NEAR(y3(0,0,2),sinc(1.005),doublethreshold);
+    EXPECT_NEAR(y3(0,0,3),sinc(2.005),doublethreshold);
+    EXPECT_NEAR(y3(0,0,4),sinc(3.005),doublethreshold);
+    EXPECT_NEAR(y3(0,0,5),sinc(3.505),doublethreshold);
+}
+
+TEST_F(bessell0,PolynomialInterpolationThreeDimension)
+{
+    double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,0.1);
+    double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,1.005);
+    double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,2.005);
+    double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,3.005);
+    double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(table3,0,0,TableLength,interval,3.505);
+
+    EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
+    EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
+    EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
+    EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
+    EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
+}
+
+TEST_F(bessell0,PolynomialInterpolationFourDimension)
+{
+    double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,0.1);
+    double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,1.005);
+    double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,2.005);
+    double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,3.005);
+    double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(table4,0,0,0,TableLength,interval,3.505);
+
+    EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
+    EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
+    EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
+    EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
+    EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
+}
+
+TEST_F(bessell0,PolynomialInterpolation)
+{
+    double y1 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,0.1);
+    double y2 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,1.005);
+    double y3 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,2.005);
+    double y4 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,3.005);
+    double y5 = ModuleBase::PolyInt::Polynomial_Interpolation(tabley,TableLength,interval,3.505);
+
+    EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
+    EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
+    EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
+    EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
+    EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
+}
+
+TEST_F(bessell0,PolynomialInterpolationXY)
+{
+    double y1 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,0.1);
+    double y2 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,1.005);
+    double y3 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,2.005);
+    double y4 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,3.005);
+    double y5 = ModuleBase::PolyInt::Polynomial_Interpolation_xy(tablex,tabley,TableLength,3.505);
+
+    EXPECT_NEAR(y1,sinc(0.1),doublethreshold);
+    EXPECT_NEAR(y2,sinc(1.005),doublethreshold);
+    EXPECT_NEAR(y3,sinc(2.005),doublethreshold);
+    EXPECT_NEAR(y4,sinc(3.005),doublethreshold);
+    EXPECT_NEAR(y5,sinc(3.505),doublethreshold);
+}