Hi Sentaro21,
This might be a good starting point, it's a Gauss-25points-Kronrod-51points integration implementation that I ported from HP41 INTDIFFEQ Module by Jean-Marc Baillard.
Before running KRONROD, the function must be loaded as a program, similar to this:
Filename:FX
sin X→X /whatever function, it must end saving the result in X
Return
------------------------------------------------------------------------------
Filename:KRONROD
10→N<CR> /Number of subintervals
"A"?→A<CR> /Lower limit
"B"?→B<CR> /Upper limit
/Load 64 nodes and weights
{64,1}->Dim Mat A<CR>
1.231760537267155E-1->Mat A[1,1]<CR>
6.158081806783294E-2->Mat A[2,1]<CR>
6.147118987142532E-2->Mat A[3,1]<CR>
6.154448300568508E-2->Mat A[4,1]<CR>
1.222424429903100E-1->Mat A[5,1]<CR>
6.112850971705305E-2->Mat A[6,1]<CR>
1.228646926107104E-1->Mat A[7,1]<CR>
6.053945537604586E-2->Mat A[8,1]<CR>
1.837189394210489E-1->Mat A[9,1]<CR>
1.194557635357848E-1->Mat A[10,1]<CR>
5.972034032417406E-2->Mat A[11,1]<CR>
2.438668837209884E-1->Mat A[12,1]<CR>
5.868968002239421E-2->Mat A[13,1]<CR>
3.030895389311078E-1->Mat A[14,1]<CR>
1.148582591457116E-1->Mat A[15,1]<CR>
5.743711636156783E-2->Mat A[16,1]<CR>
3.611723058093878E-1->Mat A[17,1]<CR>
5.595081122041232E-2->Mat A[18,1]<CR>
4.178853821930377E-1->Mat A[19,1]<CR>
1.085196244742637E-1->Mat A[20,1]<CR>
5.425112988854549E-2->Mat A[21,1]<CR>
4.730027314457150E-1->Mat A[22,1]<CR>
5.236288580640748E-2->Mat A[23,1]<CR>
5.263252843347192E-1->Mat A[24,1]<CR>
1.005359490670506E-1->Mat A[25,1]<CR>
5.027767908071567E-2->Mat A[26,1]<CR>
5.776629302412230E-1->Mat A[27,1]<CR>
4.798253713883671E-2->Mat A[28,1]<CR>
6.268100990103174E-1->Mat A[29,1]<CR>
9.102826198296365E-2->Mat A[30,1]<CR>
4.550291304992179E-2->Mat A[31,1]<CR>
6.735663684734684E-1->Mat A[32,1]<CR>
4.287284502017005E-2->Mat A[33,1]<CR>
7.177664068130844E-1->Mat A[34,1]<CR>
8.014070033500102E-2->Mat A[35,1]<CR>
4.008382550403238E-2->Mat A[36,1]<CR>
7.592592630373576E-1->Mat A[37,1]<CR>
3.711627148341554E-2->Mat A[38,1]<CR>
7.978737979985001E-1->Mat A[39,1]<CR>
6.803833381235692E-2->Mat A[40,1]<CR>
3.400213027432934E-2->Mat A[41,1]<CR>
8.334426287608340E-1->Mat A[42,1]<CR>
3.079230016738749E-2->Mat A[43,1]<CR>
8.658470652932756E-1->Mat A[44,1]<CR>
5.490469597583519E-2->Mat A[45,1]<CR>
2.747531758785174E-2->Mat A[46,1]<CR>
8.949919978782754E-1->Mat A[47,1]<CR>
2.400994560695322E-2->Mat A[48,1]<CR>
9.207471152817016E-1->Mat A[49,1]<CR>
4.093915670130631E-2->Mat A[50,1]<CR>
2.043537114588284E-2->Mat A[51,1]<CR>
9.429745712289743E-1->Mat A[52,1]<CR>
1.684781770912830E-2->Mat A[53,1]<CR>
9.616149864258425E-1->Mat A[54,1]<CR>
2.635498661503214E-2->Mat A[55,1]<CR>
1.323622919557168E-2->Mat A[56,1]<CR>
9.766639214595175E-1->Mat A[57,1]<CR>
9.473973386174152E-3->Mat A[58,1]<CR>
9.880357945340772E-1->Mat A[59,1]<CR>
1.139379850102629E-2->Mat A[60,1]<CR>
5.561932135356714E-3->Mat A[61,1]<CR>
9.95556969790498E-1->Mat A[62,1]<CR>
1.987383892330316E-3->Mat A[63,1]<CR>
9.992621049926098E-1->Mat A[64,1]<CR>
0→D<CR>
0→E<CR>
A→H<CR>
(B-A)/(2*N)→G<CR>
Do<CR>
G+H→H<CR>
64→I<CR>
1→F<CR>
Do<CR>
G*Mat A[I]→J<CR>
H-J→X<CR>
Prog "FX"<CR>
X→K<CR>
H+J→X<CR>
Prog "FX"<CR>
X+K→J<CR>
I-1→I<CR>
Mat A[I]*J+E→E<CR>
-F→F<CR>
If F>=0 Then <CR>
I-1→I<CR>
Mat A[I]*J+D→D<CR>
IfEnd<CR>
I-1→I<CR>
LpWhile I>2<CR>
H→X<CR>
Prog "FX"<CR>
X→J<CR>
Mat A[I]*J+E→E<CR>
J*Mat A1+D→D<CR >
G+H→H<CR>
N-1→N<CR>
LpWhile N>0<CR>
G*D→D<CR>
G*E→E<CR>
D<Disp> /Shows Gauss21-Integral result
E<Disp> /Shows Kronrod51-Integral result
------------------------------------------------------------------------------
Hope it helps.
Best regards.