This routine will approximate the first derivative of a function at a point in one of two ways. A quick four point polynomial estimate may be made using a step size that is provided by the user, or a more precise adaptive procedure may be use which automatically searches for the optimal step size. The code originates PPC ROM Module [1] released by the international group People Programming Computers on December 1981.
Instructions
In the adaptive routine setting a flag allows the user to view convergence of the optimal step size. The adaptive procedure delivers an error estimate along with the derivative, but 6.5 decimal digits is all that can be trusted in any case.
Both versions of FD sample on only one side of the evaluation point so that FD may be used for onesided derivatives. Discontinuities or singularities may be avoided by selecting an increment with the appropriate sign. The routine may also be used to compute partial derivatives.
 First select a SIZE which must be a minimum of 018 depending on the number of actual additional registers required by the function.
 The function must be programmed as a subroutine in program memory with a global label name of six of less characters. The input to this subroutine, namely X, will come from a data register pointed to by R11. The output, namely ƒ(x), should be left in the Xregister. ƒ(x) should not modify registers R11 through R17, flag F09 or flag F10.
 Store the global function name from step 2 in register R10.
 Store the register to hold X in the ƒ(x) subroutine in R11. This register may normally be any register other than any of R10 though R17.
 If the calculation of ƒ'(x) is desired, store the initial Xvalue “a” in the register named in step 4.
 Store an initial step size in register R12.
 Set or reset F09 to select for receptively the quick approximation, or the derivative procedure. Note that the initial step size should be large enough to the algorithm to iterate at least three times before terminating.
 If the adaptive procedure was selected in step 7 then the convergence may be viewed as an option (set F10 to do so).
 XEQ “FD”. An estimate of the first derivative will be left in the Xregister. In addition, if the adaptive procedure was used an error estimate will appear in the Yregister.
Examples

Approximate ƒ'(x) where ƒ(x) = 3.x^{3} – 4.x^{2} + 5.x + 6
:
Establish SIZE018., and place the subroutine listed below in program memory.Store the function name “EX1” in register R10.01 LBL "EX1" 02 RCL IND 11 03 ENTER^ 04 ENTER^ 05 ENTER^ 06 3 07 * 08 4 09  10 * 11 5 12 + 13 * 14 6 16 RTN
Set flag F09 to select the quick approximation method for this example.
Store 18 in R11 (pointer to register which holds the X value)
Store .1 in R12 (step size).
Store 2 in R18 (X value)XEQ “FD” the gives the answer 25, this is exact because the polynomial formula is exact for polynomials of degree three or less. 
Find the gradient of ƒ(x,y) = (x + ln(y))^{2} at the point P(2,1)
:
The gradient is simply a vector whose components are the partial derivatives of ƒ(x,y), so this example will require the calculation of two partial derivatives.
:
Create ƒ(x,y) in program memory with x in register R20 and yin register R21.The following auxiliary program will be used to setup and execute FD.01 LBL "FXY2" 02 RCL 20 03 RCL 21 04 LN 05 + 06 X^2 07 RTN
After executing GRAD the gradient will be in memory with dƒ/dx in R22 (4.000000000) and dƒ/dy in R23 (3.999999379).08 LBL "GRAD" 09 "FXY2" 10 ASTO 10 11 20 12 STO 11 13 .1 14 STO 12 15 2 16 STO 20 17 1 18 STO 21 19 CF 09 20 SF 10 21 XEQ "FD" 22 STO 22 23 ISG 11 24 ABS 25 XEQ "FD" 26 STO 23 27 RTN
For more examples refer to [1].
Method
The quick 4point polynomial approximation uses the following formula:
ƒ'(x) = [11.ƒ(x) + 18.ƒ(x+h) – 9.ƒ(x+2.h) + 2.ƒ(x+3.h)] / 6.h
where h is the current step size.
The above formula may be found in the Handbook of Mathematical Formulas, table 25.2. The formula samples the user’s function on only one size of the point where the derivative is being evaluated; on the right if h>0 and on the left if h<0. Thus a discontinuity or singularity may be avoided by selecting an increment with the appropriate sign. The formula will be exact for polynomials of degree 3 or less, and the average of estimates will be exact for polynomials of degree 4 or less. For most functions, however, the error of the polynomial is of the order of h to the 4^{th} power.
The adaptive routine works by evaluating the polynomial formula for successively smaller h values. The sequence of estimates D_{i} should be monotonic and the sequence  D_{i+1} – D_{i}  should be monotonic and decreasing. If either monotonicity condition is violated, numerical truncation is causing error in D_{i+1} and D_{i} is delivered to the user as the final output.  D_{i+1} – D_{i}  is placed in the stack as an estimate of the error. Note however that 6 1/2 digits is the most you can depend on regardless of the error estimate. The final h value is left in R12 so that any subsequent derivative evaluations at nearby points may be made by the quick polynomial formula using the nearby optimum increment value h.
Listing
Available through the repository

Requires
 XFunctions module on the HP41cv

Available as
 source code
 raw code for the V41 emulator
 bar code for the HP Wand (HP82153A)

Note
 For line by line comments refer to [1,2,3].
01 LBL "FD" 02 FS? 09 ; quick polynomial? 03 GTO 03 04 17 ; store display mode in R17 (PPC SDfunction) 05 SIGN ; for details refer to [3] 06 RDN 07 RCL d 08 STO M 09 >"\00\00" 10 X<> M 11 "*" 12 X<> M 13 STO N 14 ASTO IND L 15 RDN 16 SCI 1 ; change display mode 17 2 E3 ; R14 is the loop counter 18 STO 14 19 LBL 00 ; both part of init and part of main loop 20 RCL 12 ; see [1] 21 .7 22 * 23 RDN 24 STO 12 25 XEQ 08 26 ENTER^ 27 X<> 16 28  29 ENTER^ 30 FS? 10 31 VIEW X 32 X<> 15 33 ISG 14 34 GTO 00 35 LASTX ; main loop cont'd 36 RDN ; see [1] 37 X=0? 38 GTO 02 39 / ; test monotonicity 40 E ; see [1] 41 X<>Y 42 X<0? 43 GTO 01 44 X<Y? 45 GTO 00 46 LBL 01 ; setup stack contents 47 X<> L ; see [1] 48 LBL 02 49 R^ 50 .7 51 ST/ 12 ; restore the previous h value 52 CLX ; without lifting the stack 53 17 ; recall display mode from R17 "PPC RDfunction" 54 SIGN ; for details refer to [2] 55 ARCL IND L 56 RDN 57 RCL d 58 STO N 59 >"**" 60 X<> O 61 STO N 62 >"*****" 63 X<> N 64 STO d 65 RDN 66 CLA 67 RTN ; return 68 LBL 03 ; 4point polynomial estimate 69 . 70 STO 13 71 XEQ IND 10 72 11 73 XEQ 09 74 18 75 XEQ 09 76 9 77 XEQ 09 78 ST+ X 79 RCL 13 80  81 RCL 12 82 3 83 * 84 ST IND 11 85 ST+ X 86 / 87 ENTER^ 88 RTN 89 LBL 09 90 * 91 ST+ 13 92 RCL 12 93 ST+ IND 11 94 GTO IND 10 95 END
References
[1]  FD – First Derivative John Kennedy, August 1981 PPC ROM User’s Manual, page 144147 

[2]  SD – Store Display Mode Leigh Borkman and Keith Jarret, August 1981 PPC ROM User’s Manual, page 400 

[3]  RD – Recall Display Mode Leigh Borkman and Keith Jarret, August 1981 PPC ROM User’s Manual, page 374375 