diff --git a/.gitignore b/.gitignore
index 3e3461a..942356c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
 *.asv
+Romberg.m
diff --git a/Legendre.m b/Legendre.m
new file mode 100644
index 0000000..da9daba
--- /dev/null
+++ b/Legendre.m
@@ -0,0 +1,39 @@
+function res=Legendre(x_value_list,iterate_times)
+% \param first_x_value: start position
+% \param iterate_times: max N num of PN
+% 
+% \return res: the root value of Legendre
+% ------------------------------------------------------------%
+% P1(X)=1
+% P2(X)=X
+% ......
+% pN+2(X)= (2*n + 3) / (n + 2) * x * PN+1(x) - (n + 1) / (n + 2) * PN(x)
+% ------------------------------------------------------------%
+% example:
+% >> Legendre([0.93 -0.93 0.66 -0.66 0.23 -0.23],7);
+%
+% return 
+% the answer is
+% 0.932470
+% -0.932470
+% 0.661209
+% -0.661209
+% 0.238619
+% -0.238619
+% ------------------------------------------------------------%
+
+syms x;
+coefficient_list=sym(zeros(0,iterate_times));
+coefficient_list(1)=1;
+coefficient_list(2)=x;
+
+for index=1:iterate_times-2
+    n=index-1;
+    coefficient_list(index+2)= simplify(...
+        (2 * n + 3) / (n + 2) * x * coefficient_list(index+1) ...
+        - (n + 1) / (n + 2) * coefficient_list(index)   ...
+        );
+end
+disp(transpose(coefficient_list));
+res=newton_iteration(x_value_list,coefficient_list(iterate_times));
+end
\ No newline at end of file
diff --git a/README.md b/README.md
index a7765e4..2780a39 100644
--- a/README.md
+++ b/README.md
@@ -10,13 +10,18 @@
 |列主元高斯消去法|                       | `GE_equ.m` |
 |改进平方根法|对称矩阵的三角分解| `LDL_equ.m` |
 |Household方法|QR分解| `Household_QR.m` <br> `Household_equ.m` |
+|龙贝格积分| | `Romberg_calculus.m` |
+
 
 ## 迭代解法
 |算法名称   |  文件名  |
 |:--------:|:--------:|
 |共轭梯度法 |`CG_equ.m`|
+|牛顿迭代法 |`newton_iteration.m`|
+|四级RK法   |`Runge_Kutta.m`|
+
 
 ## 插值
 |  算法名称   |             |     文件名      |
 |:----------:|:-----------:|:---------------:|
-|三次样条插值 |第一种边界条件 <br> 第二种边界条件 <br> 第三种边界条件|`Spline1_inter.m` <br> `Spline2_inter.m` <br> `Spline3_inter.m` |
\ No newline at end of file
+|三次样条插值 |第一种边界条件 <br> 第二种边界条件 <br> 第三种边界条件|`Spline1_inter.m` <br> `Spline2_inter.m` <br> `Spline3_inter.m` |
diff --git a/Romberg_calculus.m b/Romberg_calculus.m
new file mode 100644
index 0000000..74c5700
--- /dev/null
+++ b/Romberg_calculus.m
@@ -0,0 +1,43 @@
+function res=Romberg_calculus(fx,a,b)
+% \param fx: syms expression
+% \param a:  left  bound
+% \param b:  right bound
+% \return res: the fx calculus
+% can change precision at the following
+eps=1e-6;
+buf=zeros(15);   % size
+% ------------------------------------------------------------%
+% example:
+% fx=@(x)x^2*exp(x);
+% a=0;
+% b=0.8;
+% Romberg_calculus(fx,a,b);
+%
+% return
+% the answer is
+%     0.5697
+%     0.3172
+%     0.3146
+%     0.3146
+%     0.3146
+% ------------------------------------------------------------%
+i=1;
+buf(1,1)=(b-a)/2*(fx(b)+fx(a));
+while (1)
+    i=i+1;
+    buf(1,i)=1/2*buf(1,i-1);
+    for j=1:2^(i-2)
+        buf(1,i)=buf(1,i)+(b-a)/(2^(i-1))*fx(a+(2*j-1)*(b-a)/2^(i-1)); % get T_1_i+1 value
+    end
+    
+    for m=2:i
+        k=i+1-m;
+        buf(m,k)=(4^(m-1)*buf(m-1,k+1)-buf(m-1,k))/(4^(m-1)-1);
+    end
+    if (abs((buf(i,1)-buf(i-1,1)))<eps) 
+        break;
+    end
+end
+fprintf('%.6f\n',buf(:,1));
+res=buf(i,1);
+end
\ No newline at end of file
diff --git a/Runge_Kutta.m b/Runge_Kutta.m
new file mode 100644
index 0000000..c9e4660
--- /dev/null
+++ b/Runge_Kutta.m
@@ -0,0 +1,52 @@
+function res=Runge_Kutta(fx,a,b,alpha)
+% \param fx: syms expression
+% \param a:  left  bound
+% \param b:  right bound
+% \param alpha:  y0 value
+% \return res: y_value list
+%
+% can change iterate times at the following
+N=20;
+x_list=zeros(1,15);   % size
+y_list=zeros(1,15);   % size
+% ------------------------------------------------------------%
+% example:
+% syms x y;
+% fx=y;
+% a=0;
+% b=1;
+% alpha=1;
+% Runge_Kutta(fx,a,b,alpha);
+%
+% return
+% the answer is
+%          0         0.1000    0.2000    0.3000    0.4000    0.5000  
+%          1.0000    1.1052    1.2214    1.3499    1.4918    1.6487
+
+%          0.6000    0.7000    0.8000    0.9000    1.0000         0         
+%          1.8221    2.0138    2.2255    2.4596    2.7183         0         
+        
+% ------------------------------------------------------------%
+
+h=(b-a)/N;
+syms x y;
+
+k1=h*fx;
+k2=simplify(h*subs(subs(fx,x,x+h/2),y,y+k1/2));
+k3=simplify(h*subs(subs(fx,x,x+h/2),y,y+k2/2));
+k4=simplify(h*subs(subs(fx,x,x+h),y,y+k3));
+
+x_list(1)=a;
+y_list(1)=alpha;
+
+coefficient=simplify(y+1/6*(k1+2*k2+2*k3+k4));
+disp([k1,k2,k3,k4,coefficient]);
+for i=2:N+1
+    x_list(i)=x_list(1)+(i-1)*h;
+    y_list(i)=subs(subs(coefficient,x,x_list(i-1)),y,y_list(i-1));
+end
+disp([x_list;y_list]);
+res= y_list;
+end
+
+    
\ No newline at end of file
diff --git a/newton_iteration.m b/newton_iteration.m
new file mode 100644
index 0000000..e27acd0
--- /dev/null
+++ b/newton_iteration.m
@@ -0,0 +1,55 @@
+function res=newton_iteration(x_value_list,fx)
+% \param x_value_list: list of start position
+% \param fx: syms expression
+% \return res: the root value of fx
+%
+% can change precision and max_iteration_times and min_diff_value at the following
+precision=1e-6;
+max_iteration_times=20;
+min_diff_value=1e-4;
+% ------------------------------------------------------------%
+% example:
+% syms x;
+% fx=x^2-2*x*exp(-x)+exp(-2*x);
+% newton_iteration([0.5 0.6 0.7],fx);
+%
+% return
+% the answer is
+%     0.5671
+%     0.5672
+%     0.5672
+% ------------------------------------------------------------%
+
+syms x;
+diff_fx=diff(fx);
+fx_next=x-fx/diff_fx;
+res=zeros(1,size(x_value_list,2));
+for i = 1:size(x_value_list,2)
+    first_x_value=x_value_list(i);
+    x_before=double(first_x_value);
+    x_next=double(subs(fx_next,x,x_before));
+    
+    iterate_times=1;
+    while (abs(x_next-x_before)>precision) && (iterate_times<=max_iteration_times)
+        % if f'x < min_diff_value,
+        % return
+        if abs(subs(diff_fx,x,x_before))<min_diff_value
+            fprintf('\n diff less than min_diff_value \n min diff value is %06f\n', ...
+                double(abs(subs(diff_fx,x,x_before))));
+            break;
+        end
+        
+        % print each iterate x_next value
+        fprintf('\nthe current x_value is %f,\nthe current iterate_times is %d\n', ...
+            x_next,iterate_times);
+        
+        % next iterate
+        x_before=x_next;
+        x_next=double(subs(fx_next,x,x_before));
+        iterate_times=iterate_times+1;
+    end
+    res(i)=x_next;
+end
+fprintf('\nthe answer is\n');
+fprintf('\n%f\n',res);
+end