# hp41programs

Double Exp Integration

# Numerical Integration:  Double Exponential Formulas for the HP41

Overview

1°)  Tanh-Sinh Quadrature  §-1+1 f(x) dx

1-a)  With M-Code Routines SINH  COSH  TANH  X#Y??  ( 2 programs )
1-a)  Without M-Code Routine  ( 2 programs )

2°)  Tanh-Sinh Quadrature  §ab f(x) dx

2-a)  With M-Code Routines SINH  COSH  TANH  X#Y??  ( 2 programs )
2-b)  Without M-Code Routine  ( 2 programs )

3°)  Sinh-Sinh Quadrature  §-infinity+infinity f(x) dx

3-a)  With M-Code Routines SINH  COSH  X#Y??
3-b)  Without M-Code Routine

4°)  Exp-Sinh Quadrature   §0+infinity f(x) dx

4-a)  With M-Code Routines SINH  COSH  X#Y??
4-b)  Without M-Code Routine

-The 1st versions of these programs employ M-Code routines which are listed in "A Few M-Code Routines for the HP-41"
-The 2nd versions avoid M-Code and calculate hyperbolic functions with the built-in  E^X and  E^X-1

1°)  Tanh-Sinh Quadrature  I = §-1+1 f(x) dx

a)  With M-Code Routines SINH  COSH  TANH  X#Y??

Formula:   For a given "small" number h ,    I ~  SUMk=...-2,-1,0,1,2,...   wk f(xk)

where   xk = Tanh ( Sinh ( k.h ) )                                             ( The original formulas involve a factor PI / 2 that is omitted here )
and    wk = h Cosh ( k.h ) / Cosh2 ( Sinh ( k.h )

Data Registers:           •  R00 = Function name               ( Register R00 is to be initialized before executing "TSI" )

R01 = I = R07                          R03 = k
R02 = h , h/2 , h/4 , ....              R04 to R06: temp
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

 01 LBL "TSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 STO 06  07 CLX  08 STO 03  09 XEQ IND 00  10 RCL 02  11 *  12 STO 01 13 GTO 01  14 LBL 00  15 2  16 ST/ 01  17 ST/ 02  18 STO 06  19 SIGN  20 CHS  21 STO 03          22 LBL 01  23 RCL 06  24 ST+ 03 25 RCL 02  26 RCL 03  27 *  28 SINH  29 STO 05  30 TANH  31 STO 04  32 XEQ IND 00  33 X<> 04  34 CHS  35 XEQ IND 00  36 RCL 04 37 +  38 RCL 02  39 ST* Y  40 RCL 03  41 *  42 COSH  43 *  44 RCL 05          45 COSH  46 X^2  47 /  48 RCL 01 49 +  50 STO 01  51 LASTX  52 X#Y?  53 GTO 01          54 VIEW X  55 X<> 07  56 X#Y??  57 GTO 00  58 RCL 07  59 END

( 87 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

Example:     I = §-1+1 sqrt(1-x^2) dx

-The subroutine may be

 01  LBL "T"  02  X^2  03  SIGN  04  LASTX  05  -  06  SQRT  07  RTN  08  END

"T"  ASTO 00   and if you choose  h = 0.2

0.2  XEQ "TSI"  >>>>   1.570796325                          The HP-41 displays the successive approximations
1.570796326                                                                                                       ---Execution time = 109s---

-The exact value is PI / 2

Notes:

-Thogh it will work in many cases, the criterion lines 52-53 is obviously not very good:
-If you choose h = 1 and if  f(-1) + f(1) = 0, the program would stop with a wrong result

-A better - though not perfect - criterion could be:

1-Add    DSE 08  GTO 02     after line 53
2-Add    3  STO 08  LBL 02  after line 22

3 or another positive integer...

-If the function f has a singularity for x = +/-1 , 1 - x  will not calculated with a good precision and the program may fail.

>>>  In the following variant, the HP41 uses  1 - x = [ exp ( - sinh x ) } / cosh ( sinh x )  if x > 0

Data Registers:           •  R00 = Function name              ( Register R00 is to be initialized before executing "TSI" )

R01 = I = R07                          R03 = k
R02 = h , h/2 , h/4 , ....              R04 to R09: temp
Flags: /
Subroutine:   A program that takes x in X-register ,
1 - x in Y-register ,
1 + x in Z-register and returns f(x) in X-register

 01 LBL "TSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 ENTER  07 STO 06  08 0  09 STO 03  10 XEQ IND 00  11 RCL 02  12 *  13 STO 01  14 GTO 01  15 LBL 00 16 2  17 ST/ 01  18 ST/ 02  19 STO 06  20 SIGN  21 CHS  22 STO 03  23 LBL 01  24 RCL 06  25 ST+ 03  26 RCL 02          27 RCL 03  28 *  29 SINH  30 STO 05 31 TANH  32 ENTER  33 STO 04  34 1  35 +  36 STO 08  37 RCL 05  38 CHS  39 E^X  40 RCL 05          41 COSH  42 STO 05  43 /  44 STO 09  45 R^ 46 XEQ IND 00  47 X<> 04  48 CHS  49 RCL 08  50 RCL 09  51 X<> Z  52 XEQ IND 00  53 RCL 04  54 +  55 RCL 02  56 ST* Y  57 RCL 03  58 *  59 COSH  60 * 61 RCL 05  62 X^2  63 /  64 RCL 01  65 +  66 STO 01  67 LASTX           68 X#Y?  69 GTO 01  70 VIEW X  71 X<> 07  72 X#Y??  73 GTO 00  74 RCL 07  75 END

( 104 bytes / SIZE 010 )

 STACK INPUT OUTPUT X h > 0 I

Example1:     I = §-1+1 1 / sqrt(1-x^2) dx

-Load this subroutine in main memory

 01  LBL "T"  02  RDN  03  *  04  SQRT  05  1/X  06  RTN  07  END

"T"  ASTO 00   and if you choose  h = 0.2

0.2  XEQ "TSI"  >>>>   3.141592655                          The HP-41 displays the successive approximations
3.141592654                                                                                                        ---Execution time = 2m40s---

-The exact value is PI

Example2:     I = §-1+1 1 / sqrt(1-x^4) dx

-Since  1 - x^4 = ( 1 - x ) ( 1 + x ) ( 1 + x^2 )   we can load the following subroutine in main memory

 01  LBL "T"  02  X^2  03  SIGN  04  LASTX  05  +  06  *  07  *  08  SQRT  09  1/X  10  RTN  11  END

"T"  ASTO 00   and if you choose  h = 0.2

0.2  XEQ "TSI"  >>>>   2.622057557                          The HP-41 displays the successive approximations
2.622057557                                                                                                       ---Execution time = 2m52s---

Note:

-The criterion of lines 68-69 may be improved like with the 1st version:

1-Add    DSE 08  GTO 02     after line 69
2-Add    3  STO 08  LBL 02  after line 23

b)  Without M-Code Routine

-The same results may be obtained without M-Code:

1st Version:

 01 LBL "TSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 STO 06  07 CLX  08 STO 03  09 XEQ IND 00  10 RCL 02          11 *  12 STO 01  13 GTO 01  14 LBL 00  15 2  16 ST/ 01 17 ST/ 02  18 STO 06  19 SIGN  20 CHS  21 STO 03  22 LBL 01  23 RCL 06  24 ST+ 03  25 RCL 02          26 RCL 03  27 *  28 E^X-1  29 LASTX  30 CHS  31 E^X-1  32 - 33 STO 05  34 E^X-1  35 STO Y  36 2  37 +  38 /  39 STO 04  40 XEQ IND 00  41 X<> 04  42 CHS  43 XEQ IND 00  44 RCL 04  45 +  46 RCL 02  47 ST* Y  48 RCL 03 49 *  50 E^X  51 ENTER  52 1/X  53 +  54 ST+ X  55 *  56 RCL 05          57 2  58 /  59 E^X  60 ENTER  61 1/X  62 +  63 X^2  64 / 65 RCL 01  66 +  67 STO 01          68 LASTX  69 X#Y?  70 GTO 01  71 VIEW X  72 X<> 07  73 RND  74 X<>Y  75 RND  76 X#Y?  77 GTO 00  78 RCL 07  79 END

( 104 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

-Here, the precision is controlled by the display format

-Same remark as above about the criterion of lines 69-70

2nd Version:

 01 LBL "TSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 ENTER  07 STO 06  08 0  09 STO 03  10 XEQ IND 00  11 RCL 02  12 *  13 STO 01  14 GTO 01  15 LBL 00  16 2  17 ST/ 01  18 ST/ 02  19 STO 06  20 SIGN 21 CHS  22 STO 03  23 LBL 01  24 RCL 06  25 ST+ 03  26 RCL 02  27 RCL 03  28 *  29 E^X-1  30 LASTX  31 CHS  32 E^X-1  33 -  34 STO 05          35 2  36 /  37 X<> 05  38 E^X-1  39 STO Y  40 2 41 +  42 /  43 STO 04  44 RCL 05  45 CHS  46 E^X  47 RCL X  48 RCL 05  49 E^X  50 +  51 2  52 /  53 STO 05          54 /  55 STO 09  56 1  57 R^  58 +  59 STO 08  60 X<> Z 61 XEQ IND 00  62 X<> 04  63 CHS  64 RCL 08  65 RCL 09  66 X<> Z  67 XEQ IND 00  68 RCL 04  69 +  70 RCL 02  71 ST* Y  72 RCL 03  73 *  74 E^X  75 ENTER  76 1/X  77 +  78 *  79 RCL 05  80 X^2 81 ST+ X  82 /  83 RCL 01          84 +  85 STO 01  86 LASTX  87 X#Y?  88 GTO 01  89 VIEW X  90 X<> 07  91 RND  92 X<>Y  93 RND  94 X#Y?  95 GTO 00  96 RCL 07  97 END

( 104 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h > 0 I

-Here, the precision is also controlled by the display format

-Same remark as above about the criterion of lines 87-88

2°)  Tanh-Sinh Quadrature  I = §ab f(x) dx

a)  With M-Code Routines SINH  COSH  TANH  X#Y??

-We can use the programs above with the change of arguments   x = t ( b - a ) / 2 + ( a + b ) / 2
-But your HP-41 can do the job for you.

Data Registers:           •  R00 = Function name             ( Register R00 is to be initialized before executing "TSI2" )

R01 = I = R07                          R03 = k                      R08 = ( a + b ) / 2
R02 = h , h/2 , h/4 , ....              R04 to R06: temp       R09 = ( b - a ) / 2
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

 01 LBL "TSI2"  02 STO 02  03 RDN  04 STO 09          05 X<>Y  06 ST- 09  07 +  08 2  09 ST/ 09  10 /  11 STO 08  12  E99  13 STO 07  14 SIGN  15 STO 06  16 CLX 17 STO 03  18 RCL 08  19 XEQ IND 00  20 RCL 02  21 *  22 RCL 09  23 *  24 STO 01  25 GTO 01  26 LBL 00  27 2  28 ST/ 01  29 ST/ 02  30 STO 06  31 SIGN  32 CHS 33 STO 03  34 LBL 01  35 RCL 06  36 ST+ 03  37 RCL 02  38 RCL 03  39 *  40 SINH  41 STO 05  42 TANH  43 RCL 09  44 *  45 STO 04  46 RCL 08  47 +  48 XEQ IND 00 49 X<> 04  50 CHS  51 RCL 08  52 +  53 XEQ IND 00  54 RCL 04  55 +  56 RCL 02  57 ST* Y  58 RCL 03  59 *  60 COSH  61 *  62 RCL 09  63 *  64 RCL 05 65 COSH  66 X^2  67 /  68 RCL 01          69 +  70 STO 01  71 LASTX  72 X#Y?  73 GTO 01  74 VIEW X  75 X<> 07  76 X#Y??  77 GTO 00  78 RCL 07  79 END

( 110 bytes / SIZE 010 )

 STACK INPUT OUTPUT Z a / Y b / X h I

Example:   I = §-23 exp ( -x^2 ) dx

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  RTN  06  END

"T"  ASTO 00   and if you choose  h = 0.2

-2  ENTER^
3  ENTER^
0.2  XEQ "TSI2"  >>>>   1.768288737                          The HP-41 displays the successive approximations
1.768288738                                                                                                       ---Execution time = 123s---

Notes:

-Like in the 1st paragraph, the termination criterion of lines 72-73 may be improved:

1-Add    DSE 10  GTO 02     after line 73
2-Add    3  STO 10  LBL 02  after line 34

-If the function f has a singularity for x = a  and/or  x = b ,  x - a and/or  b - x  will not calculated with a good precision and the program may fail.

>>>  In the following variant, the HP41 uses  1 - x = [ exp ( - sinh x ) } / cosh ( sinh x )  if x > 0

Data Registers:           •  R00 = Function name             ( Register R00 is to be initialized before executing "TSI2" )

R01 = I = R07                          R03 = k                     R10 = ( b - a ) / 2
R02 = h , h/2 , h/4 , ....              R04 to R09: temp      R11 = ( a + b ) / 2
Flags: /
Subroutine:   A program that takes x in X-register ,
b - x in Y-register ,
x - a in Z-register and returns f(x) in X-register

 01 LBL "TSI2"  02 STO 02  03 RDN  04 STO 10  05 X<>Y  06 ST- 10  07 +  08 2  09 ST/ 10  10 /  11 STO 11          12  E99  13 STO 07  14 SIGN  15 STO 06  16 CLX  17 STO 03  18 RCL 10  19 ENTER  20 R^ 21 XEQ IND 00  22 RCL 02  23 *  24 RCL 10  25 *  26 STO 01  27 GTO 01  28 LBL 00  29 2  30 ST/ 01  31 ST/ 02  32 STO 06  33 SIGN  34 CHS  35 STO 03  36 LBL 01  37 RCL 06  38 ST+ 03  39 RCL 02  40 RCL 03 41 *  42 SINH  43 STO 05  44 TANH  45 RCL 10  46 *  47 STO 04  48 RCL 11  49 +  50 RCL 04  51 RCL 10  52 +  53 STO 08          54 RCL 05  55 CHS  56 E^X  57 RCL 05  58 COSH  59 STO 05  60 / 61 RCL 10  62 *  63 STO 09  64 R^  65 XEQ IND 00  66 X<> 04  67 CHS  68 RCL 11  69 +  70 RCL 08  71 RCL 09  72 X<> Z  73 XEQ IND 00  74 RCL 04  75 +  76 RCL 02  77 ST* Y  78 RCL 03  79 *  80 COSH 81 *  82 RCL 10  83 *  84 RCL 05  85 X^2  86 /  87 RCL 01  88 +  89 STO 01          90 LASTX  91 X#Y?  92 GTO 01  93 VIEW X  94 X<> 07  95 X#Y??  96 GTO 00  97 RCL 07  98 END

( 130 bytes / SIZE 012 )

 STACK INPUT OUTPUT Z a / Y b / X h > 0 I

Example:     I = §28 1 / sqrt(-x^4+20.x^3-150.x^2+500.x-544) dx = §28 1 / sqrt((8-x).(x-2).(x^2-10.x+34)) dx

-Load this subroutine in main memory

 01  LBL "T"  02  STO 13  03  10  04  -  05  RCL 13  06  *  07  34  08  +  09  *  10  *  11  SQRT  12  1/X  13  RTN  14  END

"T"  ASTO 00   and if you choose  h = 0.2

2    ENTER^
8    ENTER^
0.2  XEQ "TSI2"  >>>>   0.874019185                          The HP-41 displays the successive approximations
0.874019185                                                                                                        ---Execution time = 3m23s---

Note:

-The criterion of lines 68-69 may be improved like with the 1st version:

1-Add    DSE 12  GTO 02     after line 92
2-Add    3  STO 12  LBL 02  after line 36

b)  Without M-Code Routine

-The same results may be obtained without M-Code:

1st Version:

 01 LBL "TSI2"  02 STO 02  03 RDN  04 STO 09  05 X<>Y  06 ST- 09  07 +  08 2  09 ST/ 09  10 /  11 STO 08  12  E99  13 STO 07  14 SIGN  15 STO 06  16 CLX  17 STO 03  18 RCL 08  19 XEQ IND 00  20 RCL 02 21 *  22 RCL 09  23 *  24 STO 01  25 GTO 01  26 LBL 00  27 2  28 ST/ 01  29 ST/ 02  30 STO 06  31 SIGN  32 CHS  33 STO 03  34 LBL 01  35 RCL 06          36 ST+ 03  37 RCL 02  38 RCL 03  39 *  40 E^X-1 41 LASTX  42 CHS  43 E^X-1  44 -  45 STO 05  46 E^X-1  47 STO Y  48 2  49 +  50 /  51 RCL 09  52 *  53 STO 04  54 RCL 08  55 +  56 XEQ IND 00  57 X<> 04  58 CHS  59 RCL 08  60 + 61 XEQ IND 00  62 RCL 04  63 +  64 RCL 02  65 ST* Y  66 RCL 03      67 *  68 E^X  69 ENTER  70 1/X  71 +  72 ST+ X  73 *  74 RCL 09  75 *  76 RCL 05  77 2  78 /  79 E^X  80 ENTER 81 1/X  82 +  83 X^2  84 /  85 RCL 01  86 +  87 STO 01          88 LASTX  89 X#Y?  90 GTO 01  91 VIEW X  92 X<> 07  93 RND  94 X<>Y  95 RND  96 X#Y?  97 GTO 00  98 RCL 07  99 END

( 127 bytes / SIZE 010 )

 STACK INPUT OUTPUT Z a / Y b / X h I

-Here, the precision is controlled by the display format

-Same remark as above about the criterion of lines 89-90

2nd Version:

 01 LBL "TSI2"  02 STO 02  03 RDN  04 STO 10  05 X<>Y  06 ST- 10  07 +  08 2  09 ST/ 10  10 /  11 STO 11  12  E99  13 STO 07  14 SIGN  15 STO 06  16 CLX  17 STO 03  18 RCL 10  19 ENTER  20 R^  21 XEQ IND 00  22 RCL 02  23 *  24 RCL 10 25 *  26 STO 01  27 GTO 01  28 LBL 00  29 2  30 ST/ 01  31 ST/ 02  32 STO 06          33 SIGN  34 CHS  35 STO 03  36 LBL 01  37 RCL 06  38 ST+ 03  39 RCL 02  40 RCL 03  41 *  42 E^X-1  43 LASTX  44 CHS  45 E^X-1  46 -  47 STO 05  48 2 49 /  50 X<> 05  51 E^X-1  52 STO Y  53 2  54 +  55 /  56 RCL 10          57 *  58 STO 04  59 RCL 11  60 +  61 RCL 05  62 CHS  63 E^X  64 RCL X  65 RCL 05  66 E^X  67 +  68 2  69 /  70 STO 05  71 /  72 RCL 10 73 *  74 STO 09  75 RCL 04  76 RCL 10  77 +  78 STO 08  79 X<> Z  80 XEQ IND 00  81 X<> 04  82 CHS  83 RCL 11  84 +  85 RCL 08  86 RCL 09  87 X<> Z  88 XEQ IND 00  89 RCL 04  90 +  91 RCL 02  92 ST* Y  93 RCL 03  94 *  95 E^X  96 ENTER 97 1/X  98 +  99 * 100 RCL 10 101 * 102 RCL 05         103 X^2 104 ST+ X 105 / 106 RCL 01 107 + 108 STO 01 109 LASTX 110 X#Y? 111 GTO 01 112 VIEW X 113 X<> 07 114 RND 115 X<>Y 116 RND 117 X#Y? 118 GTO 00 119 RCL 07 120 END

( 120 bytes / SIZE 012 )

 STACK INPUT OUTPUT Z a / Y b / X h > 0 I

-Here again, the precision is controlled by the display format

-Same remark as above about the criterion of lines 110-111

3°)  Sinh-Sinh Quadrature  I = §-infinity+infinity f(x) dx

a)  With M-Code Routines SINH  COSH  X#Y??

-Here, we use the change of variable   x = Sinh ( Sinh t )         where the factor PI / 2  has again been omitted

Data Registers:           •  R00 = Function name               ( Register R00 is to be initialized before executing "SSI" )

R01 = I = R07                          R03 = k
R02 = h , h/2 , h/4 , ....              R04 to R06: temp
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

 01 LBL "SSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 STO 06  07 CLX  08 STO 03  09 XEQ IND 00  10 RCL 02  11 *  12 STO 01 13 GTO 01  14 LBL 00  15 2  16 ST/ 01  17 ST/ 02  18 STO 06          19 SIGN  20 CHS  21 STO 03  22 LBL 01  23 RCL 06  24 ST+ 03 25 RCL 02  26 RCL 03  27 *  28 SINH  29 STO 05  30 SINH  31 STO 04  32 XEQ IND 00  33 X<> 04  34 CHS  35 XEQ IND 00  36 RCL 04 37 +  38 RCL 02  39 ST* Y  40 RCL 03  41 *  42 COSH  43 *  44 RCL 05          45 COSH  46 *  47 RCL 01  48 + 49 STO 01          50 LASTX  51 X#Y?  52 GTO 01  53 VIEW X  54 X<> 07  55 X#Y??  56 GTO 00  57 RCL 07  58 END

( 86 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

Example:    Evaluate  I = §-infinity+infinity  exp ( -x^2 ) dx

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  RTN  06  END

"T"  ASTO 00   and if you choose  h = 0.2

0.2   XEQ "SSI"  >>>>     1.772418685                          The HP-41 displays the successive approximations
1.772453849
1.772453851                                                                                                        ---Execution time = 124s---

Notes:

-The exact result is sqrt(PI)
-As usual, the criterion of lines 51-52 may be improved like this:

1-Add    DSE 08  GTO 02     after line 52
2-Add    3  STO 08  LBL 02  after line 22

b)  Without M-Code Routine

-Without  SINH , COSH  and  X#Y??  it yields:

 01 LBL "SSI"  02 STO 02  03  E99  04 STO 07  05 SIGN  06 STO 06  07 CLX  08 STO 03  09 XEQ IND 00  10 RCL 02  11 *  12 STO 01  13 GTO 01  14 LBL 00  15 2  16 ST/ 01 17 ST/ 02  18 STO 06  19 SIGN  20 CHS  21 STO 03          22 LBL 01  23 RCL 06  24 ST+ 03  25 RCL 02  26 RCL 03  27 *  28 E^X-1  29 LASTX  30 CHS  31 E^X-1  32 - 33 2  34 /  35 STO 05  36 E^X-1  37 LASTX  38 CHS  39 E^X-1  40 -  41 2  42 /  43 STO 04  44 XEQ IND 00  45 X<> 04  46 CHS  47 XEQ IND 00  48 RCL 04 49 +  50 RCL 02  51 ST* Y  52 RCL 03          53 *  54 E^X  55 ENTER  56 1/X  57 +  58 *  59 RCL 05  60 E^X  61 ENTER  62 1/X  63 +  64 * 65 4  66 /  67 RCL 01  68 +  69 STO 01          70 LASTX  71 X#Y?  72 GTO 01  73 VIEW X  74 X<> 07  75 RND  76 X<>Y  77 RND  78 X#Y?  79 GTO 00  80 RCL 07  81 END

( 104 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

-Here, the precision is controlled by the display format

-Same remarks as above about the criterion of lines 71-72

4°)  Exp-Sinh Quadrature   I = §0+infinity f(x) dx

a)  With M-Code Routines SINH  COSH  X#Y??

-The change of variable is now:   x = Exp ( Sinh t )         where the factor PI / 2  is again omitted

Data Registers:           •  R00 = Function name              ( Register R00 is to be initialized before executing "ESI" )

R01 = I = R07                          R03 = k
R02 = h , h/2 , h/4 , ....              R04 to R06: temp
Flags: /
Subroutine:   A program that takes x in X-register and returns f(x) in X-register

 01 LBL "ESI"  02 STO 02  03  E99  04 STO 07  05 CLX  06 STO 03  07 SIGN  08 STO 06  09 XEQ IND 00  10 RCL 02  11 *  12 STO 01 13 GTO 01  14 LBL 00  15 2  16 ST/ 01  17 ST/ 02  18 STO 06  19 SIGN  20 CHS  21 STO 03          22 LBL 01  23 RCL 06  24 ST+ 03 25 RCL 02  26 RCL 03  27 *  28 SINH  29 E^X  30 STO 04  31 STO 05  32 XEQ IND 00  33 X<> 04  34 1/X  35 XEQ IND 00  36 RCL 05 37 /  38 RCL 04  39 RCL 05  40 *  41 +  42 RCL 02  43 ST* Y  44 RCL 03          45 *  46 COSH  47 *  48 RCL 01 49 +  50 STO 01  51 LASTX  52 X#Y?  53 GTO 01          54 VIEW X  55 X<> 07  56 X#Y??  57 GTO 00  58 RCL 07  59 END

( 85 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

Example:    Evaluate  I = §0+infinity  exp ( -x^2 ) dx

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  RTN  06  END

"T"  ASTO 00   and if you choose  h = 0.2

0.2   XEQ "ESI"  >>>>     0.886226890                          The HP-41 displays the successive approximations
0.886226926
0.886226926                                                                                                        ---Execution time = 3m48s---

Notes:

-The exact result is sqrt(PI) / 2
-As usual, the criterion of lines 52-53 may be improved like this:

1-Add    DSE 08  GTO 02     after line 53
2-Add    3  STO 08  LBL 02  after line 22

b)  Without M-Code Routine

-If you want to avoid M-Code:

 01 LBL "ESI"  02 STO 02  03  E99  04 STO 07  05 CLX  06 STO 03  07 SIGN  08 STO 06  09 XEQ IND 00  10 RCL 02  11 *  12 STO 01  13 GTO 01  14 LBL 00  15 2 16 ST/ 01  17 ST/ 02  18 STO 06  19 SIGN  20 CHS  21 STO 03  22 LBL 01  23 RCL 06          24 ST+ 03  25 RCL 02  26 RCL 03  27 *  28 E^X-1  29 LASTX  30 CHS 31 E^X-1  32 -  33 2  34 /  35 E^X  36 STO 04  37 STO 05  38 XEQ IND 00  39 X<> 04  40 1/X  41 XEQ IND 00  42 RCL 05  43 /  44 RCL 04  45 RCL 05 46 *  47 +  48 RCL 02  49 ST* Y  50 RCL 03          51 *  52 E^X  53 ENTER  54 1/X  55 +  56 *  57 2  58 /  59 RCL 01  60 + 61 STO 01  62 LASTX  63 X#Y?  64 GTO 01  65 VIEW X  66 X<> 07  67 RND  68 X<>Y  69 RND  70 X#Y?  71 GTO 00  72 RCL 07          73 END

( 96 bytes / SIZE 008 )

 STACK INPUT OUTPUT X h I

-Here, the precision is controlled by the display format

-Same remarks as above about the criterion of lines 63-64

References:

[1]  Hidetosi TAKAHASI and Masatake MORI - "Double Exponential Formulas for Numerical Integration"
[2]  David H. Bailey - "Tanh-Sinh High-Precision Quadrature"