; -------- NON LINEAR SYSTEMS ---------
;      solves f(x,y)=0 and g(x,y)=0
; ----- guillaume.tello@orange.fr -----


; procedures

#label A initial_x
#label B initial_y
#label C set_precision
#label D solve
#label A' init

; variables

#reg 00 x               ; x value
#reg 01 dx              ; x variation for next loop
#reg 02 y               ; y value
#reg 03 det             ; determinant for the system
#reg 04 dy              ; y variation for next loop
#reg 05 precision       ; precision wanted

#reg 06 temp_1          ; temporary registers
#reg 07 temp_2
#reg 08 temp_3
#reg 09 temp_4
#reg 10 temp_5
#reg 11 temp_6

#reg 12 delta           ; value 1.0001

; flags

#flag 7 error_flag

;-------------------------
; X0 A to init the X value
;-------------------------

LBL initial_x
    STO x
RTN

;-------------------------
; Y0 A to init the Y value
;-------------------------

LBL initial_y
    STO y
RTN

;-------------------------
; e C to set the precision
;-------------------------

LBL set_precision
    STO precision
RTN

;----------------------
; D to solve the system
;----------------------

LBL solve
    DO
    
            temp_1,dy,det={ SBR f(x,y) }
                        +/- STO temp_6
        
            temp_2={ SBR g(x,y) }
        
        y*{ @delta }                            ; y=y+h for df/dy and dg/dy
        
            temp_5,dx={ SBR f(x,y) }
                        INV SUM det
        
            temp_3={ SBR g(x,y) }
        
        x*{ y/{ @delta } }                      ; y back to normal and x=x+h for df/dx, dg/dx
        
            temp_4={ SBR f(x,y) }
                        INV SUM temp_5
                        SUM temp_6

            det,dy*{ SBR g(x,y) }
        
            temp_5,dx,temp_4*{ @temp_2 }
        
        x/{ @delta }                            ; x back to normal
        
            temp_6,temp_1*{ @temp_3 }
        
            det+{ @temp_5 + @temp_6 = }
        
            det*{ $1e4 }
        
        dy-{ @temp_4 } dy*{ @y }

        dx-{ @temp_1 } dx*{ @x }

        error?
        IF(on) error_flag
            BREAK                   ; if an error gets out
        ENDIF
    
        mean SUM x abs              ; computes dx/det and dy/det at the same time
        X/T  SUM y abs
        IF(x>=t)
            x%t                     ; put e=max(|dx|,|dy|) in T
        ENDIF
        RCL precision X/T pause     ; e at the display and precision wanted in T
        
    UNTIL(x<t)                      ; loop until e<precision

        INV EE RCL y X/T RCL x

RTN

;---------------
; first function
;---------------

LBL E                               ; function F(x,y) = F(@reg00, @reg02)
f(x,y):
    @x X + @y X - 9 =             ; example F(x,y) = x + y - 9
RTN


#org 227

;------------------------------------
; A' to init the delta value,
; then can be overwritten by function
;------------------------------------

LBL init
    delta={ $1.0001 }
RTN

;----------------
; second function
;----------------

LBL E'                              ; function G(x,y) = G(@reg00, @reg02)
g(x,y):
    ( 4 - @x ) X + @y X - 4 =     ; example G(x,y) = (4-x) + y - 4
RTN
