'   program DIF-fin1.bas
nomainwin
 q=202
 DIM C(q,q) :DIM YY(q):dim CC(q) :dim yy(q):dim xx(q)
 OPEN "c:\kpdat" FOR OUTPUT AS #1
 WindowWidth = 590 : WindowHeight = 500'470 270
 '                               left top length hight
 StaticText #w.1,            "A:"             ,140, 10,  45, 20
 StaticText #w.2,            "B:"             ,200, 10,  45, 20
 StaticText #w.3,            "C:"             ,260, 10,  45, 20
 StaticText #w.4,            "D:"             ,320, 10,  45, 20
 STATICTEXT #w.5,            "E:"             ,380, 10,  45, 20
 staticText #w.6,       "Xstart:"             ,380, 85,  40, 20
 StaticText #w.7,       "Ystart:"             ,380, 115 ,40, 20
 StaticText #w.8,        "dY/dX:"             ,380, 145 ,40, 20
 StaticText #w.9,        "Xend:"              ,470, 85,  30, 20
 StaticText #w.10,        "Yend:"             ,470, 115, 30, 20
 StaticText #w.11,        "dY/dX:"            ,465, 145 ,35, 20
 StaticText #w.12,  " Step Fi,Ru:"            ,100, 120, 50, 40
 STATICTEXT #w.13,  "C: ONLY RUNGE"         ,240  ,63 ,100, 15
 StaticText #w.status1, "",          60,  50 ,  370,  15'diffeq
 StaticText #w.status3, "",          10,  30,   40,  40'START HERE
 StaticText #w.status5, "",          400, 175, 200,  20'recnum
 StaticText #w.status6, "",            5, 175,  70,  20'recnum
 StaticText #w.status9,  "",         10,  70,   100,  30

 TextBox #w.eq,                      100, 80,  270,  20
 TextBox #w.n,                       120, 30,   50,  20
 TextBox #w.n1,                      180, 30,   50,  20
 TextBox #w.t,                       240, 30,   50,  20
 TextBox #w.xx,                      300, 30,   50,  20
 TextBox #w.comp1,                   360, 30,   50,  20

 TextBox #w.comp2,                    420, 80,    40,  20'xin
 TextBox #w.mol2,                     420, 110,   40,  20'yin
 TextBox #w.comp3,                    420, 140,   40,  20'dy/dx in
 TextBox #w.mol3,                     500,  80,   40,  20'xout
 TextBox #w.mol4,                     140, 120,   40,  20'runge step
 TextBox #w.comp4,                    500, 110,   40,  20'yout
 TextBox #w.comp5,                    500, 140,   40,  20'dy/dx out

 Button #w.c1, "Reset", [reset],                  UL, 30,  10, 50, 30
 Button #w.c3, "Calc-Fi",      [calc],            UL, 280,120, 60, 30
 Button #w.c3, "Calc-Ru",      [calc1],           UL, 200,120, 60, 30
 Button #w.exit, "Close",   [close.w],            UL, 500,20, 40, 30
 Button #w.add, "Add New at end",   [add.rec],    UL,  80,170, 90, 25
 Button #w.upd, "Update present",   [update.rec], UL, 180,170, 90, 25
 Button #w.prv, "Previous",         [prev.rec],   UL, 280,170, 50, 25
 Button #w.next, "Next",            [next.rec],   UL, 340,170, 50, 25

 TEXTEDITOR #w.res,                                  10,200, 550,250'HEIGHT

    open "DIFFERENTIAL EQ.: INITIAL (RUNGE) AND TWO POINT VALUE (FINITE DIFF.)" for dialog as #w
    print #w, "trapclose [close.w]"
    print #w.n, "!setfocus"

    '              **** Adress file ****
    axx=6
   'prompt " Choose Indata file(diffeq) number (1-10)";axx
    open "DIFFEQ";axx;".dat" for random as #vcl len=128
    print #w.status9, "Indata file: DIFFEQ";axx;".dat "

    print #w.status1,"D2Y/DX2 +   A*DY/DX  +    B*Y      +    C*Y^2    +       D*X     +        E      =  0"

    Field #vcl, _
    60 as EQ$, _
    5 as A11, _
    5 as B11, _
    5 as C11, _
    5 as D11, _
    5 as mol1,_
    6 as YS,_
    5 as XE,_
    5 as N1,_
    5 as YE,_
    6 as DERE,_
    5 as E11,_
    5 as XS,_
    6 as DER,

     print #w.res, "  **** Saved data **** "
    rec= lof(#vcl) / 128

    print #w.status6, "Next add at ";rec+1
    for I=1 to rec
      get #vcl,I
      print #w.res, " ";EQ$ 
    next I
    print #w.res, "  **** End of saved data **** "
    print #w.eq,"!font Times_New_Roman 10 bold"
    print #w.n,"!font Times_New_Roman 10 bold"
    print #w.n1,"!font Times_New_Roman 10 bold"
    print #w.t,"!font Times_New_Roman 10 bold"
    print #w.xx,"!font Times_New_Roman 10 bold"
    print #w.comp1,"!font Times_New_Roman 10 bold"
    print #w.comp2,"!font Times_New_Roman 10 bold"
    print #w.comp3,"!font Times_New_Roman 10 bold"
    print #w.comp4,"!font Times_New_Roman 10 bold"
    print #w.comp5,"!font Times_New_Roman 10 bold"
    print #w.mol2,"!font Times_New_Roman 10 bold"
    print #w.mol3,"!font Times_New_Roman 10 bold"
    print #w.res,"!font Times_New_Roman 12  bold"
    goto [next.rec]
    goto [loop]
    PRINT" "

[validate]'subroutine
   valid = 1
    print #w.eq, "!contents?"
    input #w.eq, EQ$
    if len(EQ$)>60 then
    notice "Error!" +chr$(13)+"Allowed string=60,  Actual=";len(EQ$)
    valid=0
    return
    end if

    print #w.n, "!contents?"
    input #w.n, A11

    print #w.n1, "!contents?"
    input #w.n1, B11

    print #w.t, "!contents?"
    input #w.t, C11

    print #w.xx, "!contents?"
    input #w.xx, D11

    print #w.comp1, "!contents?"
    input #w.comp1,E11

    print #w.comp2, "!contents?"
    input #w.comp2,XS

    print #w.comp3, "!contents?"
    input #w.comp3,DER

    print #w.comp4, "!contents?"
    input #w.comp4,YE

    print #w.comp5, "!contents?"
    input #w.comp5,DERE

    print #w.mol2, "!contents?"
    input #w.mol2,YS

    print #w.mol3, "!contents?"
    input #w.mol3,XE

    print #w.mol4, "!contents?"
    input #w.mol4,N1
  return

[display.rec]'subroutine
     print #w.eq,EQ$
     print #w.n, A11
     print #w.n1,B11
     print #w.t, C11
     print #w.xx,D11
     print #w.mol2,YS
     print #w.mol3,XE
     print #w.mol4,N1
     print #w.comp1,E11
     print #w.comp2,XS
     print #w.comp3,DER
     print #w.comp4,YE
     print #w.comp5,DERE
   return

[add.rec]
    gosub [validate]
    if valid = 0 then [loop]
    recNum = lof(#vcl) / 128 + 1 ' calc location of next record
    put #vcl, recNum
    print #w.res, "!cls"
    print #w.status5, " recNum:"; recNum;"  has been added at ";recNum
    rec= lof(#vcl) / 128
    print #w.status6, "Next add at ";rec+1
  goto [loop]

[update.rec]
    gosub [validate]
    if valid = 0 then [loop]
    put #vcl, recNum
    print #w.status5, " recNum: ";recNum;"  has been updated at ";recNum
  goto [loop]

[prev.rec]
    if recNum > 1 then
    recNum = recNum - 1
    get #vcl, recNum
    yy1=lof(#vcl) / 128
    print #w.status5, "record-Nr=";recNum;":   End of file at ";yy1
    else
    print #w.status5, " Start of file:    End of file at ";yy1
    end if
    gosub [display.rec]
  goto [loop]

[next.rec]
    if recNum < lof(#vcl) / 128 then
    recNum = recNum + 1
    get #vcl, recNum
    yy=lof(#vcl) / 128
    print #w.status5, "record-Nr=";recNum;":   End of file at ";yy
    else
    print #w.status5, "End of file:   Next add at ";recNum+1
    end if
    gosub [display.rec]
  goto [loop]

[reset]
    print #w.n,     "0";
    print #w.n1,    "0" ;
    print #w.t,     "0" ;
    print #w.xx,    "0" ;
    print #w.comp1, "0" ;
    print #w.comp2, "0" ;
    print #w.mol2,  "0" ;
    print #w.comp3, "0" ;
    print #w.mol3,  "0" ;
    print #w.comp4, "0" ;
    print #w.mol4,  "0" ;
    print #w.comp5, "0" ;
    print #w.res, "!cls" ;
    print #w.res, "**** Saved data ****"
    rec= lof(#vcl) / 128
    print #w.status6, "Next add at ";rec+1
    for I=1 to rec
      get #vcl,I
      print #w.res, " ";EQ$ 
    next I
    print #w.res, "**** End data ****"
    print #w.eq, " " ;
    print #w.res,"Start: Press next to get saved data"
    print #w.res,"Change or insert new data, followed by UPDATE or ADD"
    print #w.res,"ADD places the data at the end of the file"
    print #w.res,"Play with the buttons so you learn how to preserve your data"
    print #w.res,"Y^2-term can only be handled by Runge"
    print #w.res,"Problem with repeated calculations; clear memory by closing and start new"
    print #w.res, "Use 40 steps as a rule in the iteration. Max 200 steps can be used "
    print #w.res, "The difference eq for dy/dx is of first order"
    print #w.res, "and for d2y/dx2 of seccond order with much smaller error"
    print #w.res,"Increasing the number of steps increases accuracy"
    print #w.res,"High constant Y-value >10^10 indicates division by zero"
    print #w.res,"In case with dy/dx -->dy/dx, dy/dx-->Y and Y-->dy/dx shooting,"
    print #w.res,"100-200 steps are recommended and sufficient high accuracy will be obtained"
    print #w.res,""
  goto [loop]
 return

[loop]
    input var$
    goto [loop]

[calc]

    '*************    INPUT DATA INFORMATION    *******
    'A11=A(INPUT WINDOW)  B11=B   C11=C   D11=D  E11=E
    'XS=XSTART(INPUT WINDOW) XE=XEND   YS= YSTART
    'DER=DY/DX(START)   DERE=DY/DX(END)
    'N1=FINITE STEP
    '**************************************************
'   INDATA
    M=N1-1'  M=NUMBER OF EQUATIONS
    xin  =XS
    xend =XE
    yin  =YS
    yend =YE
    uin =DER
    uend=DERE
    pp=xend-xin
    DX=pp/N1
   'END INDATA

'  SETUP OF COEFF. MATRIX
    'd2y/dx2 + A*dy/dx + B*y + [C*y^2] +D*x + E = 0
    'Difference equations j=2
    'Used eq for d2y/dx2: d2y/dx2=  (y3  -2*y2+y1 ) /dx^2
    'Used eq for dy/dx:     dy/dx=A*(y3/2     -y1/2)/dx
    'Used eq for y:          y   =        B*y2
    'y1*(1-A*dx/2) + y2*(-2+B*dx^2) + y3*(1+A*dx/2)= -(D*dx*j+E)*dx^2

                    '  C(i,j)
                    '  j=    1   2   3   4   5          6
                    '
                    'i=1    |A   1            | |Y1| |CC-YS|
                    '  2    |1   A   1        | |Y2| |CC   |
                    '  3    |    1   A   1    |*|Y3|=|CC   |
                    '  4    |        1   A   1| |Y4| |CC   |
                    '  5    |            1   A| |Y5| |CC-YE|
    n1=M
    for i=1 to n1
      xx(i)=i/(n1+1)*pp+XS
      CC(i)=(-1)*(E11+D11*(xx(i)))*DX^2
      C(i,n1+1)=CC(i)             'CONSTANT TERM will change with new x-terms
    next i
      C(1,n1+1) =C(1,n1+1)-YS
     C(n1,n1+1)=C(n1,n1+1)-YE

    for i=1 to n1-1
      C(i,i+1)=1+A11*DX/2          'UPPER START WITH y2 no change
    next i

    for i=1 to n1
       C(i,i)=  B11*DX^2-2         'DIAGONAL will change with new y-terms
    next i'

    for i=2 to n1
      C(i,i-1)=1-A11*DX/2          'LOWER START WITH y1 no change
    next i
    'END OF SETUP
    'FIRST EQUATION  AT  n1=1 giving DY/DX
     CONFIRM "Are all boundaries known?"+chr$(13)+"If answer is no use at least 100 steps "+Chr$(13)+" "+Chr$(13)+ "yes                           no";answer$ 
     IF answer$="yes" THEN
     gosub [exact]
     else
     gosub [appr]
     end if

    print #w.res,EQ$
    print #w.res,"XSTART=";XS;"   XEND=";XE;"   YSTART=";YS;"   YEND=";YE
    PRINT #w.res,"DY/DX(START)=";DER;"     DY/DX(END)=";DERE
    print #w.res,"  X               Y                 DY/DX(appr)"
    print #w.res,"     ";XS;"            ";YS;"                ";DER;"             Boundary"
    gosub [gauss]
    print #w.res,"     "; pp/(n1+1)*(n1+1)+XS;"           ";YE;"              ";DERE;"       Boundary"
    GOTO [loop]

[calc1]
        '*************    INPUT DATA INFORMATION    *******
        'A11=A(INPUT WINDOW)  B11=B   C11=C   D11=D  E11=E
        'XS=XSTART(INPUT WINDOW) XE=XEND   YS= YSTART
        'DY/DX(START)=DER   DY/DX(END)=DERE
        'N1= STEP IN RUNGE
        '**************************************************
        PRINT #w.res,""
        print #w.res,EQ$ 
        print #w.res,"XSTART= "; XS;"  XEND= ";XE
        PRINT #w.res,"YSTART ="; YS;"  YEND= ";YE
        PRINT #w.res,"USTART =";DER;"  UEND= ";DERE
        N=N1
        PRint #w.res,"RUNGE KUTTA STEPS= ";N

           XIN=XS
          XEND=XE
           YIN=YS
           YEND=YE
           UIN=DER
           UEND=DERE
         print #w.res,"  X               Y                     DY/DX"

         DX=(XEND-XIN)/N1/10
         X=XIN:  Y=YIN:  U=UIN: Z=ZIN
       PRINT"     ";X;"       ";Y;"          ";U
         PRINT #w.res,"     ";USING("###.###",X+.0001);"       ";Y;"          ";U
         FOR KA=1 TO N1
          X=XIN:  Y=YIN:  U=UIN: Z=ZIN

          FOR I=1 TO 10
             GOSUB [RUNGE]
          NEXT I
        PRINT"   ";X;"       ";Y;"        ";U
          PRINT #w.res,"   ";USING("###.###",X+.0001);"       ";Y;"        ";U
           YIN=Y
           UIN=U
           XIN=X
           NEXT KA
            print #w.res,"Runge Kutta should be used whenever possible very accurate"
            print #w.res,"Finite diff. has low accuracy for some combinations in this program"
        GOTO [loop]

[close.w]
    close #w :CLOSE #1:close #vcl
    print" press ENTER to close"
    input r$
    print"     *** end ***"
    end

[gauss]
     N=M+1
  FOR J=1 TO M
       FOR I=J TO M
         IF C(I,J)<>0 then
         FOR K=1 TO N
         X=C(J,K)
         C(J,K)=C(I,K)
         C(I,K)=X
       NEXT K
       end if
     NEXT I

     C(J,J)=C(J,J)+10^(-12)
     Y=1/C(J,J)
     FOR K=1 TO N
       C(J,K)=Y*C(J,K)
     NEXT K

     FOR I=1 TO M
       IF I=J GOTO 41
       Y=C(I,J)*(-1)
         FOR K=1 TO N
           C(I,K)=C(I,K)+Y*C(J,K)
         NEXT K
41   NEXT I 'JUMP
  NEXT J
    print""
     FOR J=1 TO M
     if J<M-1 then YY(J)=(-1)*((3)*C(J,N)-4*C(J+1,N)+1*C(J+2,N))/2/DX
     if J>2 then yy(J)=(1)*(C(J-2,N)-4*C(J-1,N)+3*C(J,N))/2/DX
     'print J;"  der= "; YY(J);"     ";yy(J)
     NEXT M
  FOR J=1 TO M
  IF J=1             THEN PRINT #w.res,using("##.###",J/(n1+1)*pp+XS);"     ";using("#####.####",C(J,N));"      ";using("#####.####",YY(J))
  IF J=2             THEN PRINT #w.res,using("##.###",J/(n1+1)*pp+XS);"     ";using("#####.####",C(J,N));"      ";using("#####.####",YY(J))
  IF J>=3 AND J< M-1 THEN PRINT #w.res,using("##.###",J/(n1+1)*pp+XS);"     ";using("#####.####",C(J,N));"      ";using("#####.####",YY(J))
  IF J=M-1           THEN PRINT #w.res,using("##.###",J/(n1+1)*pp+XS);"     ";using("#####.####",C(J,N));"      ";using("#####.####",yy(J))
  IF J=M             THEN PRINT #w.res,using("##.###",J/(n1+1)*pp+XS);"     ";using("#####.####",C(J,N));"      ";using("#####.####",yy(J))
   REM
 NEXT J

    for I=1 to N
      for K=1 to N
        C(I,J)=0
      next K
    next I
  return

[EQ]       '  SUBROUTINE EQUATIONS
        '*************    INPUT DATA INFORMATION    *******
        'A11=A(INPUT WINDOW)  B11=B   C11=C   D11=D  E11=E
        'XS=XSTART(INPUT WINDOW) XE=XEND   YS= YSTART
        'DY/DX(START)=DER   DY/DX(END)=DERE
        'N1=STEP IN RUNGE
        'd2y/dx2 + A*dy/dx + B*y + C*y^2 +D*x + E = 0
        '**************************************************
        REM  WRITTEN AS TWO FIRST ORDER: A=U:  B=(A11*U+Y*B11+C11*Y^2+D11*X+E11)*(-1)

        'A=DY B=DU  D=DZ
        '**************************************
         A=U
         B=(A11*U+Y*B11+C11*Y^2+D11*X+E11)*(-1)
        '**************************************
     RETURN

[RUNGE]      ' SUBROUTINE RUNGE
        rem  VARIABLES IN SUBROUTINE: INDEPENTENT: X  DEPENDENT:  Y, U AND Z
        X=X: Y=Y: U=U: Z=Z
        GOSUB [EQ]
        K1=A*DX:    L1=B*DX:    F1=D*DX

        X=X+DX/2: Y=Y+K1/2: U=U+L1/2: Z=Z+F1/2
        GOSUB [EQ]
        K2=A*DX:    L2=B*DX:    F2=D*DX
        X=X-DX/2: Y=Y-K1/2: U=U-L1/2 :Z=Z-F1/2

        X=X+DX/2: Y=Y+K2/2: U=U+L2/2 :Z=Z+F2/2
        GOSUB [EQ]
        K3=A*DX:    L3=B*DX:    F3=D*DX
        X=X-DX/2: Y=Y-K2/2: U=U-L2/2 :Z=Z-F2/2

        X=X+DX:   Y=Y+K3:  U=U+L3    :Z=Z+F3
        GOSUB [EQ]
        K4=A*DX:    L4=B*DX:    F4=D*DX
        X=X-DX:   Y=Y-K3:  U=U-L3    :Z=Z-F3

        X=X+DX
        Y=Y+K1/6+K2/3+K3/3+K4/6:   U=U+L1/6+L2/3+L3/3+L4/6
        Z=Z+F1/6+F2/3+F3/3+F4/6
      RETURN

[appr]  'SUBROUTINE
        'FIRST EQUATION  AT  n1=1 giving DY/DX
         print #w.res,""
         print #w.res,""
         CONFIRM "Shooting from START"+chr$(13)+"Shooting y-->y gives best fit"+chr$(13)+"Boundary  at START must be set else 0"+Chr$(13)+" "+Chr$(13)+ "DY/DX=yes              Y=no";answer$
         IF answer$="yes" THEN
         C(1,n1+1)=DER
         C(1,3)=0.5/DX*(-1)
         C(1,2)=2.0/DX*( 1)
         C(1,1)=1.5/DX*(-1)
         PRINt #w.res,"Approximative solution"
         PRINt #w.res,"Shooting from dy/dx at start";"  steps=";N1
         else
         PRINt #w.res,"Shooting from Ystart ";"  steps= ";N1
         END IF

        'LAST EQUATION  AT  n1 giving DY/DX
         CONFIRM"Shooting towards END"+chr$(13)+"Shooting y-->y gives best fit"+chr$(13)+"Boundary at END must be set else 0"+Chr$(13)+" "+Chr$(13)+    "DY/DX=yes               Y=no";answer$
         IF answer$="yes" THEN
         C(n1,n1+1)=DERE
         C(n1,n1)  =1.5/DX*( 1)
         C(n1,n1-1)=2.0/DX*(-1)
         C(n1,n1-2)=0.5/DX*( 1)
         PRINt #w.res,"Shooting towards dy/dx at end";"  steps=";N1
         else
         PRINt #w.res,"Shooting towards Yend ";"  steps= ";N1
         ee=ee+1
         END IF
      return

[exact] 'SUBROUTINE
       ' FIRST EQUATION  AT  n1=1 giving DY/DX
         print #w.res,""
         print #w.res,""
         CONFIRM "Shooting from START"+chr$(13)+"Shooting y-->y gives best fit"+chr$(13)+"Boundary  at START must be set else 0"+Chr$(13)+" "+Chr$(13)+ "DY/DX=yes              Y=no";answer$
         IF answer$="yes" THEN
         C(1,n1+1)=DER +1.5/DX*YS
         C(1,2)=0.5/DX*(-1)
         C(1,1)=2.0/DX*( 1)
         'C(1,1)=1.5/DX*(-1)
         PRINt #w.res,"Exact solution"
         PRINt #w.res,"Shooting from dy/dx at start";"  steps=";N1
         else
         PRINt #w.res,"Shooting from Ystart ";"  steps= ";N1
         END IF

        'LAST EQUATION  AT  n1 giving DY/DX
         CONFIRM"Shooting towards END"+chr$(13)+"Shooting y-->y gives best fit"+chr$(13)+"Boundary at END must be set else 0"+Chr$(13)+" "+Chr$(13)+    "DY/DX=yes               Y=no";answer$
         IF answer$="yes" THEN
         C(n1,n1+1)=DERE -1.5/DX*YE     '-2.02 used during test
         'C(n1,n1)  =1.5/DX*( 1)
         C(n1,n1-0)=2.0/DX*(-1)
         C(n1,n1-1)=0.5/DX*( 1)
         PRINt #w.res,"Shooting towards dy/dx at end";"  steps=";N1
         else
         PRINt #w.res,"Shooting towards Yend ";"  steps= ";N1
         END if
       return