The Museum of HP Calculators

# Bessel Functions for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)
2°) Bessel Functions Jn(x) of integer order ( n >= 0 )  ( improved )
3°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( non-integer order )
4°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( integer order )
5°) Continued Fractions for Jn(x) & Yn(x)  ( new )
6°) An Asymptotic Expansion for  Jn(x) and Yn(x)
7°) An Asymptotic Expansion for  Kn(x)
8°) Integrals of  Jn(x) &  In(x)  ( new )
9°) Integral of  Jn(x)  ( integer order( improved )
10°) Spherical Bessel Functions  ( new )

a) Spherical Bessel Functions of the first kind
b) Spherical Bessel Functions of the second ( and first ) kind

1°) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)

Formulae:       Jn(x) = (x/2)n [ 1/Gam(n+1)  +  (-x2/4)1/ (1! Gam(n+2) )  + .... + (-x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....

In(x) = (x/2)n [ 1/Gam(n+1)  +  (x2/4)1/ (1! Gam(n+2) )  + .... + (x2/4)k/ (k! Gam(n+k+1) ) + ....  ]      n # -1 ; -2 ; -3 ; ....

Data Registers:  /
Flag:  F01
Subroutine:  "GAM"  ( Gamma Function )

-Synthetic registers N & O may be replaced by any unused data registers.

01  LBL "JNX"
02  2
03  /
04  STO N
05  X<>Y
06  STO O
07  1
08  ENTER^
09  STO M
10  ENTER^
11  LBL 01
12  X<> Z
13  FC? 01
14  CHS
15  RCL N
16  X^2
17  *
18  RCL M
19  ST/ Y
20  RCL O
21  +
22  /
23  ISG M
24  CLX
25  ST+ Y
26  X<> Z
27  X#Y?
28  GTO 01
29  RCL N
30  RCL O
31  Y^X
32  *
33  X<> O
34  1
35  +
36  XEQ "GAM"
37  RCL O
38  X<>Y
39  /
40  CLA
41  END

( 70 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y n f(x) X x f(x)

where  f(x) = Jn(x)  if flag F01 is clear  and  f(x) = In(x)  if flag F01 is set.

Example:    Calculate   J0.7(1.9)  and  I0.7(1.9)

CF 01
0.7   ENTER^
1.9   XEQ "JNX"   yields   J0.7(1.9) = 0.584978102   (  in 12 seconds )

SF 01
0.7   ENTER^
1.9   XEQ "JNX"   yields   I0.7(1.9) = 1.727630598   (  in 12 seconds )

Notes:

-If  n  is a positive integer,  J-n(x) = (-1)n Jn(x)   and    I-n(x) = In(x)
J0(0) = 1 but this program doesn't work if x = n = 0 ( DATA ERROR line 31 )
if you add  X=Y?  X#0?  GTO 00  SIGN  X<>Y  LBL 00  after line 30 , it will work in this case too.

Jn(x) is not accurately computed for large arguments, "JNX1" is better:

2°) Bessel Functions Jn(x) of integer order ( n >= 0 )

-The following program is only a modification of a program given by Keith Jarett in "HP-41 Extended Functions made easy"
-It uses the relations:   J0(x) + 2 J2(x) + 2 J4(x) + 2 J6(x) + ...... = 1  and   Jn-1(x) + Jn+1(x) = (2n/x) Jn(x)

Data Registers:  R00 and R02 to R05: temp   R01 = x
Flags: /
Subroutines: /

01  LBL "JNX1"
02  X#0?                           Lines 02 to 08 are useful to yield  J0(0) = 1  and  Jn(0) = 0  if  n # 0
03  GTO 00
04  X#Y?
05  RTN
06  SIGN
07  RTN
08  LBL 00
09  STO 01
10  ABS
11  5
12  +
13  X<>Y
14  STO 00
15  X<Y?
16  X<>Y
17  INT
18  ST+ X
19  STO 03
20  ST- 00
21  CLST
22  STO 04
23  STO 05
24  SIGN
25  LBL 01
26  STO Z
27  RCL 03
28  ST+ X
29  *
30  RCL 01
31  /
32  X<>Y                     The sums in register R04-R05 may exceed  9.999999999 1099
33  -                              So, in order to avoid an OUT OF RANGE,
34  ISG 00                    add    ENTER^   ABS   RCL 06   X>Y?   GTO 01   ST/ 02   ST/ Z   ST/ T   ST/ 04   ST/ 05   LBL 01   R^  R^   after line 33
35  STO 02                   and add   E40   STO 06   after line 23
36  ST+ 04
37  RCL 05                   Thus you'll get, for instance:    J0(1000) = 0.02478668624  ( in 24mn26s )
38  X<> 04
39  STO 05
40  RDN
41  DSE 03
42  GTO 01
43  RCL 05
44  ST+ X
45  X<>Y
46  -
47  RCL 02
48  X<>Y
49  /
50  END

( 70 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Y n / X x Jn(x)

Examples:

2    ENTER^
10   XEQ "JNX1"   produces   J2(10) =  0.254630314  (  in 17 seconds )

3    ENTER^
100    R/S   >>>>    J3(100)  =  7.628420178 10 -2  ( in 2mn01s )

-Unlike "JNX" , "JNX1" yields accurate results for large x-values

-If you replace lines 45 thru 49 by  RCL Y   -   ST/ 02   ST/ Z   /   RCL 02
you'll get  J1(x)  in Z-register , J0(x)  in Y-register and  Jn(x)  in X-register.

3°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( non-integer order )

Formulae:     Yn(x) = ( Jn(x) cos(n(pi)) - J-n(x) ) / sin(n(pi))      ;      Kn(x) = (pi/2) ( I-n(x) - In(x) ) / sin(n(pi))      n # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....

Data Registers:   R00 = x ; R01 = n ; R02 = temp
Flag: F01      If  F01 is clear,  "YNX" calculates  Yn(x)
If  F01 is set,     "YNX" ---------   Kn(x)
Subroutine: "JNX"

01  LBL "YNX"
02  DEG
03  STO 00
04  X<>Y
05  STO 01
06  CHS
07  X<>Y
08  XEQ "JNX"
09  STO 02
10  RCL 01
11  RCL 00
12  XEQ "JNX"
13  PI
14  R-D
15  RCL 01
16  *
17  STO Z
18  FS? 01
19  CLX
20  COS
21  *
22  RCL 02
23  -
24  X<>Y
25  SIN
26  /
27  PI
28  2
29  /
30  CHS
31  FC? 01
32  ST/ X
33  *
34  END

( 54 bytes / SIZE 003 )

 STACK INPUTS OUTPUTS Y n / X x f(x)

with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01

Example:  Evaluate  Y1.4(3)  and   K1.4(3)

CF 01   1.4  ENTER^   3   XEQ "YNX"   >>>>   Y1.4(3) = 0.137821836    ( in 31 seconds )
SF 01    1.4  ENTER^   3          R/S          >>>>   K1.4(3) = 0.046088036      --------------

4°) Bessel & modified Bessel Functions of the second kind  Yn(x) &  Kn(x)  ( integer order )

Formulae:       with  psi(x) = the digamma function, we have ( positive n ):

Yn(x) = -(1/pi) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (x2/4)k + (2/pi) ln(x/2) Jn(x) - (1/pi) (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (-x2/4)k / (k!(n+k)!)

Kn(x) =  (1/2) (x/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (-x2/4)k - (-1)n ln(x/2) In(x) + (1/2) (-1)n (x/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)

Data Registers:   R00 = x/2 ; R01 = n ; R02 to R05: temp
Flag: F01      If  F01 is clear,  "YNX1" calculates  Yn(x)
If  F01 is set,     "YNX1" ---------   Kn(x)
Subroutine: "JNX"

01  LBL "YNX1"
02  STO 00
03  2
04  ST/ 00
05  X<> Z
06  STO 01
07  X<>Y
08  XEQ "JNX"
09  RCL 00
10  X^2
11  LN
12  *
13  1
14  FS? 01
15  CHS
16  RCL 01
17  Y^X
18  *
19  STO 02
20  RCL 01
21  STO 03
22  CLX
23  LBL 01
24  RCL 01
25  RCL 03
26  -
27  FACT
28  RCL 03
29  1
30  -
31  X<0?
32  CLST
33  FACT
34  ST/ Y
35  CLX
36  RCL 00
37  ST* X
38  FS? 01
39  CHS
40  LASTX
41  Y^X
42  *
43  +
44  DSE 03
45  GTO 01
46  RCL 00
47  RCL 01
48  Y^X
49  /
50  ST- 02
51  1.15443133
52  STO 05
53  RCL 01
54  LBL 02
55  X#0?
56  1/X
57  ST- 05
58  X<> L
59  DSE X
60  GTO 02
61  1
62  STO 03
63  STO 04
64  RCL 05
65  LBL 03
66  RCL 04
67  RCL 00
68  X^2
69  FC? 01
70  CHS
71  *
72  RCL 03
73  ST/ Y
74  RCL 01
75  +
76  /
77  STO 04
78  RCL 05
79  LASTX
80  1/X
81  -
82  RCL 03
83  1/X
84  -
85  STO 05
86  *
87  +
88  ISG 03
89  CLX
90  X#Y?
91  GTO 03
92  RCL 01
93  FACT
94  /
95  RCL 00
96  FS? 01
97  CHS
98  RCL 01
99  Y^X
100  *
101  RCL 02
102  +
103  FC? 01
104  PI
105  FS? 01
106  -2
107  /
108  END

( 151 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Y n / X x f(x)

with  f(x) =  Yn(x)  if  CF 01  ;   f(x) =  Kn(x)  if SF 01   ;   n = 0 , 1 , 2 , ........ , 69

Example:  Evaluate  Y2(3)  and   K2(3)

CF 01    2  ENTER^   3   XEQ "YNX1"   >>>>   Y2(3) = -0.160400393    ( in 23 seconds )
SF 01    2  ENTER^   3          R/S            >>>>    K2(3) =  0.061510458      --------------

Note:    -If n is a positive integer, Y-n(x) = (-1)n Yn(x)  and  K-n(x) = Kn(x)

5°) Continued Fractions for Jn(x) & Yn(x)

-Unlike   "JNX" &  "YNX"  ( with CF 01 )  and  "YNX1"  , the following program produces accurate results for large arguments.

Formulae:

With   p + i.q = -1/(2x) + i + (i/x) [ ( 0.52 - n2 )/( 2x + 2i + ( 1.52 - n2 )/( 2x + 4i + ..... ) ) ]
and      gn =  -1/(((2n + 2)/x) - 1/(((2n + 4)/x) - ..... ))

Jn(x)  =  sign(D) [ ( q2 + ( p - gn - n/x )2 ) ( 2q/(x.Pi) ) ] -1/2    where  D = the denominator of the second continued fraction
Yn(x) =  [ ( p - gn - n/x )/q ] Jn(x)

Data Registers:   R00 thru R13 are used by "CFRZ"  R01 = x ,  R14 = n
Flags: /
Subroutines: "CFR" & "CFRZ"  ( cf "Continued Fractions for the HP-41" )

01  LBL "JYNX"
02  STO 01
03  X<>Y
04  STO 14
05  "T"                       ( The same character as line 60 )
06  ASTO 00
07  CLST
08  RCL 01
09  XEQ "CFRZ"
10  RCL 01
11  ST/ Z
12  /
13  1
14  +
15  X<>Y
16  CHS
17  RCL 01
18  ST+ X
19  1/X
20  -
21  STO 09
22  X<>Y
23  STO 10
24  "U"                       ( The same character as line 75 )
25  ASTO 00
26  CLX
27  CLA
28  RCL 01
29  XEQ "CFR"
30  CHS
31  STO 08
32  RCL 09
33  +
34  RCL 14
35  RCL 01
36  /
37  -
38  STO 11
39  RCL 10
40  R-P
41  LAST X
42  ST+ X
43  PI
44  RCL 01
45  *
46  /
47  SQRT
48  X<>Y
49  /
50  RCL 04
51  SIGN
52  *
53  STO 12
54  RCL 11
55  *
56  RCL 10
57  /
58  RCL 12
59  RTN
60  LBL "T"
61  RCL 03
62  ST+ X
63  RCL 01
64  ST+ X
65  RCL 03
66  .5
67  -
68  X^2
69  RCL 14
70  X^2
71  -
72  0
73  X<>Y
74  RTN
75  LBL "U"
76  RCL 02
77  RCL 14
78  +
79  ST+ X
80  RCL 01
81  /
82  1
83  CHS
84  END

( 121 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS Y n >= 0 Yn(x) X x > 0 Jn(x)

Examples:

10   ENTER^  XEQ "JYNX"   >>>>   J10(10) =  0.207486107   X<>Y   Y10(10) =  -0.359814151    ( in 2mn42s )

3.14  ENTER^
100  XEQ "JYNX"  >>>>  J3.14(100) =  0.079535723     X<>Y    Y3.14(100) =  0.006582327    ( in 4mn35s )

-The method doesn't work if n is a negative integer.
-If n < 0 we can use the relations  Jn = J-n cos n.Pi + Y-n sin n.Pi  and  Yn = -J-n sin n.Pi + Y-n cos n.Pi
-If x < 0 the results are generally complex.

6°) An Asymptotic Expansion for  Jn(x) and Yn(x)

Formulae:     Jn(x) = (2/(pi*x))1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) )  and   Yn(x) = (2/(pi*x))1/2 ( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )

where   P ~  1 - (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)(4n2-49)/(4!(8x)4) - ......

and   Q ~  (4n2-1)/(8x) - (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......

Data Registers:  R00 = x ; R01 = n ; R02 to R06: temp
Flags: /
Subroutines: /

01  LBL "AEJY"
02  DEG
03  STO 00
04  X<>Y
05  STO 01
06  ST+ X
07  X^2
08  STO 04
09  1
10  STO 02
11  -
12  RCL 00
13  8
14  *
15  STO 05
16  ST* 05
17  /
18  STO 03
19  XEQ 01
20  STO 06
21  CLX
22  STO 02
23  SIGN
24  STO 03
25  XEQ 01
26  GTO 02
27  LBL 01
28  RCL 03
29  RCL 02
30  2
31  +
32  STO 02
33  ST+ X
34  3
35  -
36  X^2
37  RCL 04
38  -
39  *
40  RCL 04
41  RCL 02
42  ST+ X
43  DSE X
44  X^2
45  -
46  *
47  RCL 02
48  ST/ Y
49  DSE X
50  /
51  RCL 05
52  /
53  STO 03
54  +
55  X#Y?
56  GTO 01
57  RTN
58  LBL 02
59  RCL 00
60  RCL 01
61  ST+ X
62  1
63  +
64  4
65  /
66  PI
67  *
68  -
69  R-D
70  STO 04
71  X<>Y
72  P-R
73  RCL 04
74  RCL 06
75  P-R
76  ST+ T
77  RDN
78  -
79  2
80  PI
81  /
82  RCL 00
83  /
84  SQRT
85  ST* Z
86  *
87  END

( 112 bytes / SIZE 007 )

 STACK INPUTS OUTPUTS Y n Yn(x) X x Jn(x)

Example:       4   ENTER^  100   XEQ "AEJY"  >>>>   J4(100)  =  0.026105809
X<>Y                Y4(100)  = -0.075430120    ( in 10 seconds )

-The program will not work if  x is too "small" or if n is too "large"

7°) An Asymptotic Expansion for  Kn(x)

Formula:     Kn(x) ~  (pi/(2x))1/2 e-x  ( 1 + (4n2-1)/(8x) + (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3) + ......  )

Data Registers:  R00 = x ; R01 = 4n2 ; R02 , R03: temp
Flags: /
Subroutines: /

01  LBL "AEK"
02  STO 00
03  X<>Y
04  ST+ X
05  X^2
06  STO 01
07  SIGN
08  STO 02
09  STO 03
10  LBL 01
11  RCL 03
12  RCL 01
13  RCL 02
14  ST+ X
15  DSE X
16  X^2
17  -
18  *
19  RCL 00
20  RCL 02
21  *
22  8
23  *
24  /
25  STO 03
26  +
27  ISG 02
28  CLX
29  X#Y?
30  GTO 01
31  PI
32  RCL 00
33  ST+ X
34  /
35  SQRT
36  *
37  RCL 00
38  E^X
39  /
40  END

( 54 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Y n / X x Kn(x)

Example:       2   ENTER^  10   XEQ "AEK"  >>>>   K2(10)  = 0.000,021,509,817,05    ( in 15 seconds )
1.4  ENTER^  19       R/S          >>>>   K1.4(19) = 1.683198846 10-9          ( in 6 seconds )

-The program will not work if  x is too "small" or if n is too "large".
-For instance, n = 2 , x = 7 dont work.
-However, if we stop the iterations before the terms begin to increase, we find  K2(7) = 0.0005545622
-Divergence is more typical than convergence for asymptotic series...

Remark:     In(x) is accurately computed by "JNX" with flag F01 set, but if you need an asymptotic expansion for this function:

replace line 39 with *
line 36 with /
line 34 with *  and add CHS after line 24   ( example:   I1.4(19) = 15,597,340  ( in 8 seconds ) )

8°) Integrals of   Jn(x)  &  In(x)

-Here, we integrate the series expansions used by "JNX"     ( n # -1 ; -2 ; -3 ; .... )

Data Registers:  /
Flag:  F01      CF 01  for  §0x Jn(t).dt
SF 01 for  §0x In(t).dt
Subroutine:  "1/G+"  ( Reciprocal of the Gamma Function )

-Synthetic registers M , N & O may be replaced by any unused data registers.

01  LBL "ITIJ"
02  STO N
03  X<>Y
04  STO O
05  2
06  ST/ N
07  SIGN
08  STO M
09  +
10  /
11  ENTER^
12  ENTER^
13  LBL 01
14  X<> Z
15  FC? 01
16  CHS
17  RCL N
18  X^2
19  *
20  RCL M
21  ST/ Y
22  RCL O
23  +
24  ST/ Y
25  RCL M
26  +
27  1
28  -
29  ST* Y
30  2
31  +
32  /
33  ISG M
34  CLX
35  ST+ Y
36  X<> Z
37  X#Y?
38  GTO 01
39  RCL N
40  RCL O
41  Y^X
42  *
43  X<> O
44  1
45  +
46  XEQ "1/G+"     or    XEQ "1/G"    but   If you prefer   XEQ "GAM"  or  XEQ "GAM+"
47  RCL O                                                  replace line 48 by  X<>Y   /
48  *
49  CLA
50  END

( 85 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y n / X x §0x fn(t).dt

where   fn(t) = Jn(t)  if flag F01 is clear
and    fn(t) =  In(t)  if flag F01 is set

Examples:

CF 01
1.4   ENTER^
3     XEQ"ITIJ"  >>>>   §03  J1.4(x).dx =  1.049262785  ( in 14s )

SF 01
1.4   ENTER^
3     XEQ"ITIJ"  >>>>   §03  I1.4(x).dx =  2.918753200  ( in 13s )

-If flag F01 is clear, this program cannot be used for large arguments ( like "JNX" )

9°) Integral of   Jn(x)   ( integer order )

Formula:       §0x Jn(t).dt = 2 ( Jn+1(x) + Jn+3(x) + Jn+5(x) + ........ )

Data Registers:   R00 thru R05 ( or R06 ) are used by "JNX1" ( R01 = x )  R07: temp R08: partial sums and eventually,  (1/2) §0x Jn(t).dt
Flag: F07
Subroutine: "JNX1"

01  LBL "ITJ"
02  STO 01
03  X<>Y
04  STO 07
05  CLX
06  STO 08
07  DSE 07
08  LBL 01
09  SF 07
10  LBL 02
11  RCL 07
12  2
13  +
14  STO 07
15  RCL 01
16  XEQ "JNX1"              XEQ "JNX" would work too, provided x is not too "large"
17  RCL 08                      but in this case, "ITIJ" with CF 01 is much faster!
18  +
19  STO 08
20  LASTX
21  X#Y?
22  GTO 01
23  FS?C 07
24  GTO 02
25  ST+ X
26  END

( 45 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS Y n / X x §0x Jn(t).dt

Examples:

1    ENTER^
3    XEQ "ITJ"  >>>>  §03  J1(x).dx = 1.260051955  ( in 2mn01s )                 - "ITIJ" yields the same result in 14 seconds -

0    ENTER^
10      R/S         >>>>  §010  J0(x).dx = 1.067011304  ( in 5mn54s )                 - "ITIJ" yields 1.067011425 ( error ~ 10 -7 )  in 23 seconds -

50    ENTER^
30      R/S         >>>>  §030  J50(x).dx = 1.478729947 10 -8  ( in 12mn29s )

-For small arguments, "ITJ" is much slower than "ITIJ"
-There is also a risk of  OUT OF RANGE  for large x-values,
but if you modify "JNX1" as explained on the right of its listing ( cf paragraph 2°) )
you'll get for instance §0100  J50(x).dx = 1.088806747 ( in about 10 seconds with a V41 emulator in turbo mode )

10°) Spherical Bessel Functions

a) Spherical Bessel Functions of the first kind

-Spherical Bessel functions of the first kind are defined by   jn(x) = (Pi/(2x))1/2 Jn+1/2(x)   where n is an integer
-So, we can compute these functions by one of the programs above.
-However, another approach is used hereafter:

Formulae:             jn-1(x) =  jn(x) (2n+1)/x - jn+1(x)   this recurrence relation is used , starting with  jm(x) = 1 , jm+1(x) = 0  for some large m
-Then, the results are normalized by the sum     1 j02(x) + 3 j12(x) + 5 j22(x) + ....... = 1
-We could also normalize by calculating  j0(x) = (sin x)/x  but  this would be inaccurate if  sin x ~ 0

Data Registers:   R00 , R02 , R03: temp    R01 = x , R04 = Sum
Flags: /
Subroutines: /

01  LBL "SB1"
02  X#0?                           Lines 02 to 08 are only useful to yield   j0(0) = 1  and  jn(0) = 0  if  n # 0
03  GTO 00
04  X#Y?
05  RTN
06  SIGN
07  RTN
08  LBL 00
09  STO 01
10  ABS
11  5
12  +
13  X<>Y
14  STO 00
15  X<Y?
16  X<>Y
17  INT
18  ST+ X
19  STO 03
20  ST- 00
21  ST+ X
22  STO 04
23  CLST
24  SIGN
25  ST+ 04
26  LBL 01
27  RCL 03
28  ST+ X
29  1
30  +
31  RCL Y
32  *
33  RCL 01
34  /
35  R^                           The sum in register R04 may easily exceed  9.999999999 1099
36  -                              So, in order to avoid an OUT OF RANGE,
37  ISG 00                    add    ENTER^   ABS   RCL 05   X>Y?   GTO 01   ST/ 02   ST/ Z   ST/ T   X^2   ST/ 04   LBL 01   R^  R^   after line 36
38  STO 02                   and add   E40   STO 05  after line 23
39  ENTER^
40  X^2                         Thus you'll get, for instance,   j100(50) = 1.019012264 10 -22  ( in 3mn22s )
41  RCL 03
42  ST+ X
43  DSE X
44  *
45  ST+ 04
46  RDN
47  DSE 03
48  GTO 01
49  RCL 02
50  RCL 04
51  SQRT
52  /
53  END

( 74 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS Y n >= 0 / X x jn(x)

Examples:

2  PI  XEQ "SB1"  >>>>  j2(Pi) =  0.303963551   ( execution time = 14 seconds )

10  ENTER^
2      R/S        >>>>  j10(2) =  6.825300865 10 -8  ( in 17 seconds )

100  ENTER^  R/S   >>>>  j100(100) =  0.01088047702  ( in 3mn )

b) Spherical Bessel Functions of the second ( and first ) kind

-Spherical Bessel functions of the second kind are defined by    yn(x) = (Pi/(2x))1/2 Yn+1/2(x)   where n is an integer
-But the following program uses another method:

Formulae:        yn+1(x) =  yn(x) (2n+1)/x - yn-1(x)                                  ;       jn+1(x) =  jn(x) (2n+1)/x - jn-1(x)
y0(x) = -(cos x)/x  ,  y1(x) = -(cos x)/x2 - (sin x)/x          ;        j0(x) = (sin x)/x  ,   j1(x) = (sin x)/x2 - (cos x)/x

Data Registers:   R00 = k.nnn  ; R01 = x
Flag: F01    CF 01 to calculate  yn(x)
SF 01 -----------   jn(x)          Set flag F01 only if  x is not significantly smaller than  n  -  the reccurence relation is unstable otherwise.
Subroutines: /

01  LBL "SB2"
03  STO 00
04  X<>Y
05   E3
06  /
07  STO 01
08  SIGN
09  P-R
10  RCL 00
11  ST/ Z
12  /
13  FC? 01
14  CHS
15  DEG
16  FS? 01
17  X<>Y
18  LBL 01
19  ISG 01
20  FS? 30
21  GTO 02              this line may be replaced by RTN and line 33 may be deleted
22  STO Z
23  RCL 01
24  INT
25  ST+ X
26  DSE X
27  *
28  RCL 00
29  /
30  X<>Y
31  -
32  GTO 01
33  LBL 02
34  END

( 53 bytes / SIZE 002 )

 STACK INPUTS OUTPUTS Y 0 <= n < 1000 / X x # 0 fn(x)

where  fn(x) = yn(x)  if  CF 01
and   fn(x) = jn(x)   if  SF 01     ( only if  x is not significantly smaller than  n )

Examples:

a)         CF 01
2     ENTER^
3.14  XEQ "SB2"  >>>>   y2(3.14) = -0.222053752  ( in 2 seconds )

CF 01
30   ENTER^
5   XEQ "SB2"  >>>>   y5(30) = -7.760717577 1018  ( in 15 seconds )

b)         SF 01
4    ENTER^
100  XEQ "SB2"  >>>>   j4(100) = -4.179461836 10 -3  ( in 3 seconds )

SF 01
100  ENTER^   XEQ "SB2"  >>>>   j100(100) = 1.088047702 10 -2  ( in 47 seconds )

-If  n < 0 , we can use the relation    yn(x) = (-1)n+1 j -n-1(x)   which actually holds for all integers n  ( n < 0 , n = 0 , n > 0 )

References:     Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4