hp41programs

nthorderODE

Nth Order Differential Equations for the HP-41


Overview
 

 1°) Second order Differential Equations

   a) Explicit Method of Order 4
   b) Explicit Method of Order 6
   c) Implicit Method of Order 6

 2°) Third order Differential Equations

   a) Explicit Method of Order 4
   b) Explicit Method of Order 6

 3°) Nth order Differential Equations

   a) Explicit Method of Order 4     ( X-Functions Module required )
   b) Explicit Method of Order 6

-Nth order differential equations can be re-written as a system of first order differential equations, so these programs may seem redundant !
-They are, however, much easier to use, especially for the general case.
 

Recently added:   §1b §1c §2b & §3b
 
 

1°) Second Order Differential Equations
 

     a) Explicit Method of Order 4
 

-We want to solve the equation  y" = f(x,y,y')     with the initial values:   y(x0) = y0  ,   y'(x0) = y'0
 

Data Registers:           •  R00 = function name                       ( Registers R00 thru R05 are to be initialized before executing "SRK"  )

                                      •  R01 = x0
                                      •  R02 = y0        •  R04 = h/2   where h = stepsize
                                      •  R03 = y'0       •  R05 =  m  = number of steps             R06 thru R09: temp
Flags: /
Subroutine:       A program which computes  y" = f(x,y,y')  assuming  x , y , y'  are in registers X , Y , Z ( respectively )  upon entry
 
 

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

 
   ( 103 bytes / SIZE 010 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:     Let's consider the Lane-Emden equation of index 3     y" = -(2/x) y' - y3   with the initial values  y(0) = 1 , y'(0) = 0

-There is already a difficulty because x = 0 is a singular point!
-But if we use a series expansion, we find after substitution that  y = 1 + a.x2 + ....  will satisfy the LEE  if   a = -1/6    whence  y"(0) = -1/3

-So, we can load the following subroutine into program memory:

  LBL "T"      ST+ X      3              RTN             CHS                (  LBL "T" or any other global LBL , maximum 6 characters )
  X=0?          X<>Y      Y^X         LBL 00         RTN
  GTO 00      /               +              3
  RCL Z        X<>Y      CHS         1/X

-If we want to estimate y(1)  with a stepsize h = 0.1   ( whence  m = 10 )

  "T"  ASTO 00
   0   STO 01  STO 03             0.05  STO 04   ( h/2 )
   1   STO 02                            10    STO 05

  XEQ "SRK"   >>>>       x   =  1                            ( in 56 seconds )
                        RDN     y(1)  =  0.855057170
                        RDN     y'(1) = -0.252129561

-These new "initial" values are also stored in registers R01 R02 R03
-With  h/2 = 0.025  and  m = 20 we would have found  y(1) = 0.855057539  &  y'(1) = -0.252129276
 

Notes:    The solution of the Lane-Emden Equation of index 3 looks like this:

      y
    1|*                I
      |                  *
      |                              *
      |                                           *
      |                                                              *
      |-----------------------------------------------------------*x1-------- x
     0
              y(x1) = 0  for  x1 = 6.896848619    and    y'(x1) = -0.0424297576
and there is an inflexion point I with  xI = 1.495999168 , yI = 0.720621687 and  y'(xI) = -0.279913175

-The solutions of the Lane-Emden Equations of index n   ( y" + (2/x).y' + yn = 0  ,  y(0)= 1 , y'(0) = 0 )
  can be expressed by elementary functions for only 3 n-values:

   a)  n = 0   y(x) = 1 - x2/6
   b)  n = 1   y(x) = (sin x)/x
   c)  n = 5   y(x) = ( 1+x2/3) -1/2
 

     b) Explicit Method of Order 6
 

-"RKN6" employs a 7-stage 6th-order formula ( cf Runge-Kutta -> Runge-Kutta-Nystrom-General ODEs for the HP-41, example2 )
 
 

Data Registers:           •  R00 = function name                    ( Registers R00 thru R05 are to be initialized before executing "RKN6" )

                                      •  R01 = x          •  R04 = h                R06 to R12: temp
                                      •  R02 = y          •  R05 = N
                                      •  R03 = y'
Flags: /
Subroutine:  A program that takes  x , y , y'  in registers  X , Y , Z  respectively and returns  f(x,y,y')  in X-register.
 

-Line 327 is a three-byte  GTO 01
 
 

 01 LBL "RKN6"
 02 RCL 05
 03 STO 06
 04 LBL 01
 05 RCL 03
 06 RCL 02
 07 RCL 01
 08 XEQ IND 00
 09 RCL 04
 10 *
 11 STO 07
 12 3
 13 /
 14 RCL 03
 15 +
 16 RCL 07
 17 18
 18 /
 19 RCL 03
 20 3
 21 /
 22 +
 23 RCL 04
 24 *
 25 RCL 02
 26 +
 27 RCL 04
 28 3
 29 /
 30 RCL 01
 31 +
 32 XEQ IND 00
 33 RCL 04
 34 *
 35 STO 08
 36 ST+ X
 37 3
 38 /
 39 RCL 03
 40 +
 41 RCL 08
 42 ST+ X
 43 9
 44 /
 45 RCL 03
 46 ST+ X
 47 3		  48 /
 49 +
 50 RCL 04
 51 *
 52 RCL 02
 53 +
 54 RCL 04
 55 ST+ X
 56 3
 57 /
 58 RCL 01
 59 +
 60 XEQ IND 00
 61 RCL 04
 62 *
 63 STO 09
 64 CHS
 65 RCL 08
 66 4
 67 *
 68 +
 69 RCL 07
 70 +
 71 12
 72 /
 73 RCL 03
 74 +
 75 RCL 07
 76 RCL 09
 77 +
 78 36
 79 /
 80 RCL 03
 81 3
 82 /
 83 +
 84 RCL 04
 85 *
 86 RCL 02
 87 +
 88 RCL 04
 89 3
 90 /
 91 RCL 01
 92 +
 93 XEQ IND 00
 94 RCL 04		  95 *
 96 STO 10
 97 90
 98 *
 99 RCL 09
100 35
101 *
102 +
103 RCL 08
104 110
105 *
106 -
107 RCL 07        
108 25
109 *
110 +
111 48
112 /
113 RCL 03
114 +
115 RCL 10
116 270
117 *
118 RCL 09
119 35
120 *
121 +
122 RCL 08
123 330
124 *
125 -
126 RCL 07
127 125
128 *
129 +
130 288
131 /
132 RCL 03
133 5
134 *
135 6
136 /
137 +
138 RCL 04
139 *
140 RCL 02
141 +		 142 RCL 04
143 5
144 *
145 6
146 /
147 RCL 01
148 +
149 XEQ IND 00
150 RCL 04
151 *
152 STO 11
153 10
154 /
155 RCL 10
156 2
157 /
158 +
159 RCL 09
160 8
161 /
162 -
163 RCL 08
164 11
165 *
166 24
167 /
168 -
169 RCL 07
170 3
171 *
172 20
173 /
174 +
175 RCL 03
176 +
177 RCL 11
178 15
179 /
180 CHS
181 RCL 10
182 12
183 /
184 -
185 RCL 09
186 16
187 /
188 +		 189 RCL 08
190 11
191 *
192 144
193 /
194 +
195 RCL 07
196 40
197 /
198 +
199 RCL 03        
200 6
201 /
202 +
203 RCL 04
204 *
205 RCL 02
206 +
207 RCL 04
208 6
209 /
210 RCL 01
211 +
212 XEQ IND 00
213 RCL 04
214 *
215 STO 12
216 1600
217 *
218 RCL 11
219 128
220 *
221 +
222 RCL 10
223 2360
224 *
225 -
226 RCL 09
227 215
228 *
229 +
230 RCL 08
231 1980
232 *
233 +
234 RCL 07
235 783		 236 *
237 -
238 780
239 /
240 RCL 03
241 +
242 RCL 12
243 4 E3
244 *
245 RCL 11
246 64
247 *
248 +
249 RCL 10        
250 4720
251 *
252 -
253 RCL 09 
254 215
255 *
256 +
257 RCL 08
258 3960
259 *
260 +
261 RCL 07
262 2349
263 *
264 -
265 2340
266 /
267 RCL 03
268 +
269 RCL 04
270 *
271 RCL 02
272 +
273 RCL 04
274 RCL 01
275 +
276 XEQ IND 00
277 RCL 04
278 *
279 RCL 12
280 80
281 *
282 RCL 11
283 16		 284 *
285 +
286 RCL 10
287 110
288 *
289 +
290 RCL 09
291 55
292 *
293 +
294 RCL 07         
295 39
296 *
297 +
298 600
299 /
300 RCL 03 
301 +
302 RCL 04
303 ST+ 01
304 *
305 ST+ 02
306 CLX
307 RCL 07
308 +
309 13
310 *
311 RCL 09
312 RCL 10
313 +
314 55
315 *
316 +
317 RCL 11
318 RCL 12
319 +
320 32
321 *
322 +
323 200
324 /
325 ST+ 03
326 DSE 06
327 GTO 01
328 RCL 03
329 RCL 02
330 RCL 01
331 END

 
      ( 433 bytes / SIZE 013 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /          y'(x)
           Y             /          y(x)
           X             /            x

 
Example:     y" = - 2 y - 2 x y'          y(0) = 1    y'(0) = 0         Calculate  y(1)
 

-Load for example this routine in main memory:
 
 

01 LBL "T"
02 ST* Z
03 RDN
04 +
05 ST+ X
06 CHS
07 RTN

 
 "T"  ASTO 00

  0  STO 01  STO 03
  1  STO 02

-With  h = 0.1  &  N = 10

  0.1  STO 04  10  STO 05  XEQ "RKN6"   >>>>   x    =  1                     = R01
                                                                     RDN  y(1) = 0.367879455   = R02
                                                                     RDN y'(1) = -0.735758821  = R03

Notes:

-The exact solution is  y = exp(-x2)  so  y(1) = 0.367879441....
-Simply press R/S to continue with the same h & N
 

     c) Implicit Method of Order 6
 

 "IRKN6" employs the implicit 3-stage 6th-order formula described in reference [2]
 
 

Data Registers:           •  R00 = function name                    ( Registers R00 thru R05 are to be initialized before executing "IRKN6" )

                                      •  R01 = x          •  R04 = h                R06 to R13: temp
                                      •  R02 = y          •  R05 = N
                                      •  R03 = y'
Flags: /
Subroutine:  A program that takes  x , y , y'  in registers  X , Y , Z  respectively and returns  f(x,y,y')  in X-register.
 

-Lines 215-252 are three-byte  GTO 01
 
 

 01 LBL "IRKN6"
 02 RCL 05
 03 STO 06
 04 5
 05 STO 13
 06 15
 07 SQRT
 08 STO 07
 09 ST+ 13
 10 -
 11 10
 12 ST/ 13
 13 /
 14 STO 12
 15 CLX
 16 STO 08
 17 STO 09
 18 STO 10
 19 LBL 01
 20 RCL 08 
 21 25
 22 *
 23 40
 24 RCL 07
 25 12
 26 *
 27 -
 28 RCL 09
 29 *
 30 +
 31 25
 32 RCL 07
 33 6
 34 *
 35 -
 36 RCL 10
 37 *		  38 +
 39 180
 40 /
 41 RCL 03 
 42 +
 43 14
 44 RCL 07         
 45 4
 46 *
 47 -
 48 RCL 09
 49 *
 50 RCL 08 
 51 ST+ X
 52 +
 53 RCL 07
 54 5
 55 *
 56 20
 57 -
 58 RCL 10 
 59 *
 60 -
 61 180
 62 /
 63 RCL 03
 64 RCL 12
 65 *
 66 +
 67 RCL 04
 68 *
 69 RCL 02
 70 +
 71 RCL 04
 72 RCL 12
 73 *
 74 RCL 01		  75 +
 76 XEQ IND 00
 77 RCL 04
 78 *
 79 ENTER
 80 X<> 08
 81 -
 82 ABS
 83 STO 11
 84 RCL 07
 85 3
 86 *
 87 10
 88 +
 89 RCL 08
 90 *
 91 RCL 09
 92 16
 93 *
 94 +
 95 10
 96 RCL 07 
 97 3
 98 *
 99 -
100 RCL 10
101 *
102 +
103 72
104 /
105 RCL 03
106 +
107 7
108 RCL 07
109 ST+ X
110 +
111 RCL 08		 112 *
113 RCL 09
114 4
115 *
116 +
117 7
118 RCL 07 
119 ST+ X
120 -
121 RCL 10        
122 *
123 +
124 144
125 /
126 RCL 03
127 2
128 /
129 +
130 RCL 04
131 *
132 RCL 02
133 +
134 RCL 04
135 2
136 /
137 RCL 01
138 +
139 XEQ IND 00
140 RCL 04
141 *
142 ENTER
143 X<> 09
144 -
145 ABS
146 ST+ 11
147 25
148 RCL 07		 149 6
150 *
151 +
152 RCL 08 
153 *
154 RCL 07
155 12
156 *
157 40
158 +
159 RCL 09        
160 *
161 +
162 RCL 10
163 25
164 *
165 +
166 180
167 /
168 RCL 03 
169 +
170 RCL 07
171 5
172 *
173 20
174 +
175 RCL 08 
176 *
177 RCL 07
178 4
179 *
180 14
181 +
182 RCL 09
183 *
184 +
185 RCL 10		 186 ST+ X
187 +
188 180
189 /
190 RCL 03 
191 RCL 13
192 *
193 +
194 RCL 04
195 *
196 RCL 02
197 +
198 RCL 04
199 RCL 13
200 *
201 RCL 01
202 +
203 XEQ IND 00
204 RCL 04
205 *
206 ENTER
207 X<> 10
208 -
209 ABS
210 ST+ 11
211 RCL 11
212  E-9
213 VIEW 11
214 X<Y?
215 GTO 01
216 5
217 RCL 07
218 +
219 RCL 08
220 *
221 RCL 09		 222 8
223 *
224 +
225 5
226 RCL 07 
227 -
228 RCL 10
229 *
230 +
231 36
232 /
233 RCL 03         
234 +
235 RCL 04
236 ST+ 01
237 *
238 ST+ 02
239 RCL 08 
240 RCL 10 
241 +
242 5
243 *
244 RCL 09
245 8
246 *
247 +
248 18
249 /
250 ST+ 03
251 DSE 06
252 GTO 01
253 RCL 03
254 RCL 02
255 RCL 01
256 CLD
257 END

 
      ( 322 bytes / SIZE 014 )
 
 

      STACK        INPUTS      OUTPUTS
           Z            /          y'(x)
           Y             /          y(x)
           X             /            x

 
Example:     y" = - 2 y - 2 x y'          y(0) = 1    y'(0) = 0         Calculate  y(1)
 

-Load for example this routine in main memory:
 
 

01 LBL "T"
02 ST* Z
03 RDN
04 +
05 ST+ X
06 CHS
07 RTN

 
 "T"  ASTO 00

  0  STO 01  STO 03
  1  STO 02

-With  h = 0.1  &  N = 10

  0.1  STO 04  10  STO 05  XEQ "IRKN6"   >>>>   x    =   1                     = R01
                                                                      RDN  y(1) =  0.367879441   = R02
                                                                      RDN y'(1) = -0.735758882   = R03

Notes:

-The HP41 displays the successive differences between 2 approximations of the implicit formula ( line 213 )
-Line 212 ( E-9 ) may lead to an infinite loop: change this number if need be.

-The exact solution is  y = exp(-x2)  so  y(1) = 0.367879441....
-Simply press R/S to continue with the same h & N

-This formula is the best one listed on this page for second order ODEs, but of course not very fast with an HP41.
-The error for y(1) is only about 2 E-10 in the above example.
 

2°) Third Order Differential Equations
 

     a) Explicit Method of Order 4
 

-Likewise, we want to solve the equation  y"' = f(x,y,y',y")     with the initial values:   y(x0) = y0  ,   y'(x0) = y'0  ,  y"(x0) = y"0
 

Data Registers:           •  R00 = function name                         ( Registers R00 thru R06 are to be initialized before executing "TRK"  )

                                      •  R01 = x0
                                      •  R02 = y0
                                      •  R03 = y'0       •  R05 = h/2   where h = stepsize
                                      •  R04 = y"0      •  R06 =  m  = number of steps             R07 thru R12: temp
Flags: /
Subroutine:       A program which computes  y"' = f(x,y,y',y")  assuming  x , y , y' , y"  are in registers X , Y , Z , T ( respectively )  upon entry
 
 
 

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

 
     ( 143 bytes / SIZE 013 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+m.h)
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:      y"' =  2xy" - x2y' + y2   with  y(0) = 1 , y'(0) = 0 , y"(0) = -1

-We load, for instance:

   LBL "S"     ST+ X     -                         (  LBL "S" or any other global LBL , maximum 6 characters )
   X^2           ST* T      -
   ST* Z        RDN       RTN
   X<> L       X^2

-With  h =  0.1  and  m = 10

  "S"  ASTO 00
   0   STO 01  STO 03                  0.05  STO 05   ( h/2 )
   1   STO 02  CHS  STO 04         10    STO 06

  XEQ "TRK"   >>>>       x   =  1                            ( in 49 seconds )
                        RDN     y(1)  =  0.595434736
                        RDN     y'(1) = -0.776441445
                        RDN    y"(1) =  -0.791715205

-With  h/2 = 0.025  and  m = 20,  we would find  y(1) = 0.595431304  ,  y'(1) = -0.776444326  ,  y"(1) = -0.791718300
 

     b) Explicit Method of Order 6
 

-We can get a better precision with a higher order Runge-Kutta formula:
 

Data Registers:           •  R00 = function name                         ( Registers R00 thru R06 are to be initialized before executing "TRK6"  )

                                      •  R01 = x0
                                      •  R02 = y0
                                      •  R03 = y'0       •  R05 =  h  = stepsize
                                      •  R04 = y"0      •  R06 =  N = number of steps             R07 thru R25: temp
Flags: /
Subroutine:       A program which computes  y''' = f(x,y,y',y")  assuming  x , y , y' , y"  are in registers X , Y , Z , T ( respectively )  upon entry
 

-Line 433 is a three-byte GTO 01
 
 

 01 LBL "TRK6"
 02 RCL 06
 03 STO 25
 04 LBL 01
 05 RCL 04
 06 STO 08
 07 RCL 03
 08 STO 07
 09 RCL 02
 10 RCL 01
 11 XEQ IND 00
 12 RCL 05
 13 ST* 07
 14 ST* 08
 15 *
 16 STO 09
 17 3
 18 /
 19 RCL 04
 20 +
 21 STO 11
 22 RCL 08
 23 3
 24 /
 25 RCL 03
 26 +
 27 STO 10
 28 RCL 07
 29 3
 30 /
 31 RCL 02
 32 +
 33 RCL 05
 34 3
 35 /
 36 RCL 01
 37 +
 38 RCL 11
 39 RDN
 40 XEQ IND 00
 41 RCL 05
 42 ST* 10
 43 ST* 11
 44 *
 45 STO 12
 46 1.5
 47 /
 48 RCL 04
 49 +
 50 STO 14
 51 RCL 11
 52 1.5
 53 /
 54 RCL 03
 55 +
 56 STO 13
 57 RCL 10
 58 1.5
 59 /
 60 RCL 02
 61 +
 62 RCL 05
 63 1.5		  64 /
 65 RCL 01
 66 +
 67 RCL 14
 68 RDN
 69 XEQ IND 00
 70 RCL 05
 71 ST* 13
 72 ST* 14
 73 *
 74 STO 15
 75 CHS
 76 RCL 12
 77 4
 78 *
 79 +
 80 RCL 09
 81 +
 82 12
 83 /
 84 RCL 04
 85 +
 86 STO 17
 87 RCL 11
 88 4
 89 *
 90 RCL 14
 91 -
 92 RCL 08
 93 +
 94 12
 95 /
 96 RCL 03
 97 +
 98 STO 16
 99 RCL 10
100 4
101 *
102 RCL 13
103 -
104 RCL 07
105 +
106 12
107 /
108 RCL 02
109 +
110 RCL 05
111 3
112 /
113 RCL 01
114 +
115 RCL 17
116 RDN
117 XEQ IND 00
118 RCL 05
119 ST* 16
120 ST* 17
121 *
122 STO 18
123 90
124 *
125 RCL 15
126 35		 127 *
128 +
129 RCL 12
130 110
131 *
132 -
133 RCL 09
134 25
135 *
136 +
137 48
138 /
139 RCL 04 
140 +
141 STO 20        
142 RCL 17
143 90
144 *
145 RCL 14
146 35
147 *
148 +
149 RCL 11
150 110
151 *
152 -
153 RCL 08
154 25
155 *
156 +
157 48
158 /
159 RCL 03
160 +
161 STO 19
162 RCL 16
163 90
164 *
165 RCL 13
166 35
167 *
168 +
169 RCL 10
170 110
171 *
172 -
173 RCL 07
174 25
175 *
176 +
177 48
178 /
179 RCL 02
180 +
181 RCL 05
182 1.2
183 /
184 RCL 01
185 +
186 RCL 20
187 RDN
188 XEQ IND 00
189 RCL 05		 190 ST* 19
191 ST* 20
192 *
193 STO 21
194 10
195 /
196 RCL 18
197 2
198 /
199 +
200 RCL 15
201 8
202 /
203 -
204 RCL 12        
205 11
206 *
207 24
208 /
209 -
210 RCL 09
211 .15
212 *
213 +
214 RCL 04
215 +
216 STO 23
217 RCL 20
218 10
219 /
220 RCL 17
221 2
222 /
223 +
224 RCL 14
225 8
226 /
227 -
228 RCL 11
229 11
230 *
231 24
232 /
233 -
234 RCL 08
235 .15
236 *
237 +
238 RCL 03 
239 +
240 STO 22
241 RCL 19
242 10
243 /
244 RCL 16
245 2
246 /
247 +
248 RCL 13
249 8
250 /
251 -
252 RCL 10		 253 11
254 *
255 24
256 /
257 -
258 RCL 07
259 .15
260 *
261 +
262 RCL 02
263 +
264 RCL 05
265 6
266 /
267 RCL 01
268 +
269 RCL 23
270 RDN
271 XEQ IND 00
272 RCL 05
273 ST* 22
274 ST* 23
275 *
276 STO 24
277 40
278 X^2
279 *
280 RCL 21
281 128
282 *
283 +
284 RCL 18
285 2360
286 *
287 -
288 RCL 15
289 215
290 *
291 +
292 RCL 12
293 1980
294 *
295 +
296 RCL 09
297 783
298 *
299 -
300 780
301 /
302 RCL 04
303 +
304 X<> 11
305 1980
306 *
307 RCL 23
308 40
309 X^2
310 *
311 +
312 RCL 20
313 128
314 *		 315 +
316 RCL 17
317 2360
318 *
319 -
320 RCL 14
321 215
322 *
323 +
324 RCL 08
325 783
326 *
327 -
328 780
329 /
330 RCL 03        
331 +
332 ENTER
333 X<> 10
334 1980
335 *
336 RCL 22 
337 40
338 X^2
339 *
340 +
341 RCL 19
342 128
343 *
344 +
345 RCL 16
346 2360
347 *
348 -
349 RCL 13
350 215
351 *
352 +
353 RCL 07
354 783
355 *
356 -
357 780
358 /
359 RCL 02
360 +
361 RCL 01
362 RCL 05
363 +
364 RCL 11
365 RDN
366 XEQ IND 00
367 RCL 05
368 ST+ 01
369 ST* 10
370 ST* 11
371 *
372 STO 12
373 RCL 09
374 +
375 13
376 *		 377 RCL 15
378 RCL 18
379 +
380 55
381 *
382 +
383 RCL 21
384 RCL 24
385 +
386 32
387 *
388 +
389 .5
390 %
391 ST+ 04
392 RCL 08         
393 RCL 11
394 +
395 13
396 *
397 RCL 14 
398 RCL 17
399 +
400 55
401 *
402 +
403 RCL 20
404 RCL 23
405 +
406 32
407 *
408 +
409 .5
410 %
411 ST+ 03
412 RCL 07
413 RCL 10
414 +
415 13
416 *
417 RCL 13
418 RCL 16
419 +
420 55
421 *
422 +
423 RCL 19
424 RCL 22 
425 +
426 32
427 *
428 +
429 .5
430 %
431 ST+ 02
432 DSE 25
433 GTO 01
434 RCL 04
435 RCL 03
436 RCL 02
437 RCL 01
438 END

 
     ( 617 bytes / SIZE 026 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+N.h)
           Z             /      y'(x0+N.h)
           Y             /      y(x0+N.h)
           X             /        x0+N.h

-Simply press R/S to continue with the same h and m

Example:      y"' =  2xy" - x2y' + y2   with  y(0) = 1 , y'(0) = 0 , y"(0) = -1

-We load, for instance:

   LBL "S"     ST+ X     -
   X^2           ST* T      -
   ST* Z        RDN       RTN
   X<> L       X^2

-With  h =  0.1  and  N = 10

  "S"  ASTO 00
   0   STO 01  STO 03                   0.1   STO 05
   1   STO 02  CHS  STO 04         10    STO 06

  XEQ "TRK6"   >>>>       x   =  1
                          RDN     y(1)  =  0.595431073
                          RDN     y'(1) = -0.776444516
                          RDN    y"(1)  = -0.791718501

-With  h = 0.05  &  N = 20,  it yields  y(1) = 0.5954310715  ,  y'(1) = -0.7764445234  ,  y"(1) = -0.7917185205

Note:

-The exact results, rounded to 10 decimals are:

   y(1) = 0.5954310718  ,  y'(1) = -0.7764445229  ,  y"(1) = -0.7917185202
 

3°) N-th Order Differential Equations
 

     a) Explicit Method of Order 4
 

-The differential equation is now  y(n) = f(x,y,y',y",.....,y(n-1))   with the initial values:   y(x0) = y0  ,   y'(x0) = y'0  ,  ........  ,  y(n-1)(x0) = y(n-1)0
-"NRK" employs the "classical" fourth-order Runge-Kutta formula.
 

Data Registers:           •  R00 = function name          ( Registers R00 thru R03 and R09 thru Rn+9 are to be initialized before executing "NRK"  )

                                      •  R01 = n      ( order )
                                      •  R02 = h/2    where h = stepsize                       R04 thru R08 & Rn+10 thru R3n+9: temp
                                      •  R03 =  m  =  number of steps

                                      •  R09 = x0     •  R10 = y0     •  R11 = y'0    •  R12 = y"0   .....................  •  Rn+9 = y(n-1)0
Flags: /
Subroutine:   A program which computes   y(n) = f(x,y,y',y",.....,y(n-1))   assuming  x , y , y' , y" , ........ , y(n-1)  are in registers  R09 R10 R11 R12 .... Rn+9

-If you don't have an HP-41CX, replace line 24 with  ENTER^  CLX  LBL 00  STO IND Y  ISG Y  GTO 00
 
 

01  LBL "NRK" 
02  RCL 03
03  STO 08
04  RCL 01
05  9.009
06  +
07  STO 04
08  INT
09  RCL 01
10  +
11  STO 05
12  LASTX
13  +
14  STO 06
15  RCL 04
16  INT
17  RCL 05
18  1
19  ST+ Z		 20  +
21   E3
22  /
23  +
24  CLRGX 
25  10
26  LASTX
27  +
28  RCL 01
29   E6
30  /
31  +
32  STO 07
33  GTO 03
34  LBL 01
35  XEQ IND 00
36  LBL 02
37  RCL 02
38  *		 39  ST+ IND 05
40  FS? 05
41  ST+ IND 05
42  FC? 06
43  ST+ X
44  RCL IND 06
45  +
46  X<> IND 04
47  DSE 06
48  DSE 05
49  DSE 04
50  GTO 02
51  RCL 01
52  ST+ 04
53  ST+ 05
54  ST+ 06
55  RTN
56  LBL 03
57  RCL 07		 58  REGMOVE 
59  CF 05
60  SF 06
61  XEQ 01
62  SF 05
63  RCL 02
64  ST+ 09
65  XEQ 01
66  CF 06
67  XEQ 01
68  CF 05
69  RCL 02
70  ST+ 09
71  XEQ 01
72  RCL 07
73  REGSWAP 
74  RCL 04
75  RCL 05
76  3		 77  SIGN
78  LBL 04
79  CLX
80  X<> IND Y 
81  LASTX
82  /
83  ST+ IND Z 
84  DSE Y
85  DSE Z
86  GTO 04
87  DSE 08
88  GTO 03
89  RCL 12        
90  RCL 11
91  RCL 10
92  RCL 09
93  END

 
      ( 150 bytes / SIZE 3n+10 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+m.h)
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:      y(5) = y(4) - 2x.y"' + y" - y.y'   with  y(0) = 1 , y'(0) = y"'(0) = y(4)(0) = 0 , y"(0) = -1

-We load, for instance:

   LBL "W"       RCL 13        +                  -                        (  LBL "W" or any other global LBL , maximum 6 characters )
   RCL 14         *                  RCL 10        RTN
   RCL 09         -                  RCL 11
   ST+ X           RCL 12       *

-With  h =  0.1  and  m = 10

  "W"   ASTO 00                                                       and the initial values:     0  STO 09   STO 11  STO 13  STO 14
     5    STO 01        ( fifth order equation )                                                    1  STO 10  CHS  STO 12
  0.05  STO 02        ( h/2 = 0.05 )
   10    STO 03        ( m = 10 )

  XEQ "NRK"   >>>>       x   =  1                      = R09                  ( in 2mn48s )
                        RDN     y(1)  =  0.491724880   = R10
                        RDN     y'(1) = -1.041200697   = R11       and     RCL 13 =  y"'(1)  = -0.479803795
                        RDN    y"(1) =  -1.163353624  = R12                  RCL 14 = y""(1) = -0.897595630

-With  h/2 = 0.025  and  m = 20,  it yields:

     y(1) = 0.491724223 ,  y'(1) = -1.041200398 ,  y"(1) = -1.163353549 ,  y"'(1)  = -0.479804004 ,  y(4)(1) = -0.897594479
 

     b) Explicit Method of Order 6
 

Data Registers:           •  R00 = function name          ( Registers R00 to R03 and R19 to Rn+19 are to be initialized before executing "NRK6"  )

                                      •  R01 =  n      ( order )
                                      •  R02 =  h  =  stepsize                       R04 thru R14 & Rn+20 thru R8n+19: temp               R15-R16-R17-R18:  unused
                                      •  R03 =  m  =  number of steps

                                      •  R19 = x0     •  R20 = y0     •  R21 = y'0    •  R22 = y"0   .....................  •  Rn+19 = y(n-1)0
Flags: /
Subroutine:   A program which computes   y(n) = f(x,y,y',y",.....,y(n-1))   assuming  x , y , y' , ........ , y(n-1)  are in registers  R19 R20 R21 .... Rn+19
 

-Line 367 is a three-byte GTO 10
 
 

 01 LBL "NRK6"
 02 RCL 03
 03 STO 12
 04 RCL 01
 05 RCL 01
 06 RCL 01
 07 STO 13
 08 19.019
 09 +
 10 STO 04
 11 +
 12 STO 05
 13 +
 14 STO 06
 15 +
 16 STO 07
 17 +
 18 STO 08
 19 +
 20 STO 09
 21 +
 22 STO 10
 23 +
 24 STO 11
 25 LBL 10
 26 RCL 19
 27 STO 14
 28 LBL 01
 29 RCL IND 04
 30 STO IND 05
 31 RCL 02
 32 *
 33 DSE 06
 34 DSE 13
 35 STO IND 06
 36 DSE 05
 37 DSE 04
 38 GTO 01
 39 RCL 01
 40 ST+ 04
 41 ST+ 05
 42 ST+ 06
 43 STO 13
 44 XEQ IND 00
 45 RCL 02
 46 *
 47 STO IND 06
 48 LASTX
 49 3
 50 /
 51 RCL 14
 52 +
 53 STO 19		  54 LBL 02
 55 RCL IND 06
 56 3
 57 /
 58 RCL IND 05
 59 +
 60 STO IND 04
 61 RCL 02
 62 *
 63 DSE 07
 64 DSE 13
 65 STO IND 07
 66 DSE 06
 67 DSE 05
 68 DSE 04
 69 GTO 02
 70 RCL 01
 71 ST+ 04
 72 ST+ 05
 73 ST+ 06
 74 ST+ 07
 75 STO 13
 76 XEQ IND 00
 77 RCL 02
 78 *
 79 STO IND 07
 80 RCL 14
 81 RCL 02
 82 1.5
 83 /
 84 +
 85 STO 19
 86 LBL 03
 87 RCL IND 07
 88 1.5
 89 /
 90 RCL IND 05
 91 +
 92 STO IND 04
 93 RCL 02
 94 *
 95 DSE 08
 96 DSE 13
 97 STO IND 08
 98 DSE 07
 99 DSE 05
100 DSE 04
101 GTO 03
102 RCL 01
103 ST+ 04
104 ST+ 05
105 ST+ 07
106 ST+ 08		 107 STO 13
108 XEQ IND 00
109 RCL 02
110 *
111 STO IND 08
112 RCL 14
113 RCL 02
114 3
115 /
116 +
117 STO 19
118 LBL 04
119 RCL IND 07
120 4
121 *
122 RCL IND 08
123 -
124 RCL IND 06
125 +
126 12
127 /
128 RCL IND 05
129 +
130 STO IND 04
131 RCL 02
132 *
133 DSE 09
134 DSE 13
135 STO IND 09
136 DSE 08
137 DSE 07
138 DSE 06
139 DSE 05
140 DSE 04
141 GTO 04
142 RCL 01
143 ST+ 04
144 ST+ 05
145 ST+ 06
146 ST+ 07
147 ST+ 08
148 ST+ 09
149 STO 13
150 XEQ IND 00
151 RCL 02
152 *
153 STO IND 09
154 RCL 14
155 RCL 02
156 1.2
157 /
158 +
159 STO 19		 160 LBL 05
161 RCL IND 09
162 90
163 *
164 RCL IND 08
165 35
166 *
167 +
168 RCL IND 07
169 110
170 *
171 -
172 RCL IND 06
173 25
174 *
175 +
176 48
177 /
178 RCL IND 05
179 +
180 STO IND 04
181 RCL 02
182 *
183 DSE 10
184 DSE 13
185 STO IND 10
186 DSE 09
187 DSE 08
188 DSE 07
189 DSE 06
190 DSE 05
191 DSE 04
192 GTO 05
193 RCL 01
194 ST+ 04
195 ST+ 05
196 ST+ 06
197 ST+ 07
198 ST+ 08
199 ST+ 09
200 ST+ 10
201 STO 13
202 XEQ IND 00
203 RCL 02
204 *
205 STO IND 10
206 RCL 14
207 RCL 02
208 6
209 /
210 +
211 STO 19
212 LBL 06		 213 RCL IND 10
214 10
215 /
216 RCL IND 09
217 2
218 /
219 +
220 RCL IND 08
221 8
222 /
223 -
224 RCL IND 07
225 11
226 *
227 24
228 /
229 -
230 RCL IND 06
231 .15
232 *
233 +
234 RCL IND 05
235 +
236 STO IND 04
237 RCL 02
238 *
239 DSE 11
240 DSE 13
241 STO IND 11
242 DSE 10
243 DSE 09
244 DSE 08
245 DSE 07
246 DSE 06
247 DSE 05
248 DSE 04
249 GTO 06
250 RCL 01
251 ST+ 04
252 ST+ 05
253 ST+ 06
254 ST+ 07
255 ST+ 08
256 ST+ 09
257 ST+ 10
258 ST+ 11
259 STO 13
260 XEQ IND 00
261 RCL 02
262 *
263 STO IND 11
264 RCL 14
265 RCL 02		 266 +
267 STO 19
268 RCL IND 07
269 LBL 07
270 1980
271 *
272 RCL IND 11
273 40
274 X^2
275 *
276 +
277 RCL IND 10
278 128
279 *
280 +
281 RCL IND 09
282 2360
283 *
284 -
285 RCL IND 08
286 215
287 *
288 +
289 RCL IND 06
290 783
291 *
292 -
293 780
294 /
295 RCL IND 05
296 +
297 STO IND 04
298 RCL 02
299 *
300 DSE 07
301 DSE 13
302 X<> IND 07
303 DSE 11
304 DSE 10
305 DSE 09
306 DSE 08
307 DSE 06
308 DSE 05
309 DSE 04
310 GTO 07
311 RCL 01
312 ST+ 04
313 ST+ 05
314 ST+ 06
315 ST+ 07
316 ST+ 08
317 ST+ 09
318 ST+ 10		 319 ST+ 11
320 STO 13
321 XEQ IND 00
322 RCL 02
323 *
324 STO IND 07
325 LBL 08
326 RCL IND 06
327 RCL IND 07
328 +
329 13
330 *
331 RCL IND 08
332 RCL IND 09
333 +
334 55
335 *
336 +
337 RCL IND 10
338 RCL IND 11
339 +
340 32
341 *
342 +
343 .5
344 %
345 RCL IND 05
346 +
347 STO IND 04
348 DSE 11
349 DSE 10
350 DSE 09
351 DSE 08
352 DSE 07
353 DSE 06
354 DSE 05
355 DSE 04
356 GTO 08
357 RCL 01
358 ST+ 04
359 ST+ 05
360 ST+ 06
361 ST+ 07
362 ST+ 08
363 ST+ 09
364 ST+ 10
365 ST+ 11
366 DSE 12
367 GTO 10
368 RCL 22
369 RCL 21
370 RCL 20
371 RCL 19
372 END

 
     ( 611 bytes / SIZE 8n+20 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /      y"(x0+m.h)
           Z             /      y'(x0+m.h)
           Y             /      y(x0+m.h)
           X             /        x0+m.h

-Simply press R/S to continue with the same h and m

Example:      y(5) = y(4) - 2x.y"' + y" - y.y'   with  y(0) = 1 , y'(0) = y"'(0) = y(4)(0) = 0 , y"(0) = -1

-We load, for instance:

   LBL "W"       RCL 23        +                  -                        (  LBL "W" or any other global LBL , maximum 6 characters )
   RCL 24         *                  RCL 20        RTN
   RCL 19         -                  RCL 21
   ST+ X           RCL 22       *

-With  h =  0.1  and  m = 10

  "W"   ASTO 00

     5    STO 01        ( fifth order equation )                 and the initial values:     0  STO 19   STO 21  STO 23  STO 24
   0.1   STO 02        ( h = 0.1 )                                                                      1  STO 20  CHS  STO 22
   10    STO 03        ( m = 10 )

  SIZE 060 ( or larger )

  XEQ "NRK6"   >>>>       x   =  1                      = R19
                          RDN     y(1)  =  0.491724179   = R20
                          RDN     y'(1) = -1.041200379   = R21       and     RCL 23 =  y"'(1)  = -0.479804016
                          RDN    y"(1) =  -1.163353546  = R22                  RCL 24 = y""(1)  = -0.897594395

-With  h = 0.05  and  m = 20,  it yields:

     y(1) = 0.491724180 ,  y'(1) = -1.041200380 ,  y"(1) = -1.163353546 ,  y'''(1)  = -0.479804017 ,  y""(1) = -0.897594397
 
 
 

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  Sa Agam, Ya Yahaya, Sc Osuala - "An Implicit Runge-Kutta Method for General Second Order Ordinary Differential Equations"