PROGRAM laplace
! implementation of Jacobi relaxation method to solve
! Laplace's equation in rectangular geometry
DIM V(0:100,0:100)
CALL assign(V(,),nx,ny,min_change)
CALL iterate(V(,),nx,ny,min_change,iterations)
END

SUB assign(V(,),nx,ny,min_change)
    ! linear dimension of rectangular region
    LET nx = 9     ! number of interior points in x-direction
    LET ny = 9     ! number of interior points in y-direction
    LET V0 = 10              ! boundary potential of rectangle
    INPUT prompt "percentage change = ": min_change
    LET min_change = min_change/100
    ! fix potential on boundary of rectangle
    FOR x = 0 to nx + 1
        LET V(x,0) = V0
        LET V(x,ny+1) = V0
    NEXT x
    FOR y = 0 to ny + 1
        LET V(0,y) = V0
        LET V(nx+1,y) = V0
    NEXT y
    ! guess initial values of potentials of interior sites
    FOR y = 1 to ny
        FOR x = 1 to nx
            LET V(x,y) = 0.9*V0
        NEXT x
    NEXT y
    CALL show_output(V(,),nx,ny,0)
END SUB

SUB iterate(V(,),nx,ny,min_change,iterations)
    DIM Vave(100,100)
    LET iterations = 0
    DO
       ! maximum difference of iterated potential at any site
       LET change = 0
       FOR y = 1 to ny
           FOR x = 1 to nx
               ! average of potential of neighboring cells
               LET Vave(x,y) = V(x+1,y) + V(x-1,y)
               LET Vave(x,y) = Vave(x,y) + V(x,y+1) + V(x,y-1)
               LET Vave(x,y) = 0.25*Vave(x,y)
               ! compute percentage change in potential
               IF Vave(x,y) <> 0 then
                  LET dV = abs((V(x,y) - Vave(x,y))/Vave(x,y))
                  IF dV > change then LET change = dV
               END IF
           NEXT x
       NEXT y
       FOR y = 1 to ny       ! update potential at each site
           FOR x = 1 to nx
               LET V(x,y) = Vave(x,y)
           NEXT x
       NEXT y
       LET iterations = iterations + 1
       CALL show_output(V(,),nx,ny,iterations)
    LOOP until change <= min_change
END SUB

SUB show_output(V(,),nx,ny,iterations)      ! print potential   
    PRINT
    PRINT "iteration ="; iterations
    FOR y = 0 to ny + 1
        FOR x = nx + 1 to 0 step - 1
            PRINT using "###.##": V(x,y);
        NEXT x
        PRINT
    NEXT y
END SUB