'SOLVING OF BOUNDARY VALUE PROBLEM BY A FINITE DIFFERENCE METHOD
'DIFF.EQ.:D2Y/DX2+A*Y=B: A=1/4: B=8
'THE REGION 0<X<10 IS DIVIDED INTO n1+1 EQUALLY SPACED INTERVALS OF LENGTH DX
'AT POINT J THE EQ. CAN NOW BE REPRESENTED
'Y(J+1)-2*Y(J)+Y(J-1)/DX^(2)+A*Y(J)=B
'AA=A*DX^(2)-2
'IF n1= 19 THEN DX=0.5 AND AND A SET OF 19 EQ. CAN BE SET UP
'A TRIDIAGONAL MATRIX IS OBTAINED
'THE FIRST EQ IS [AA  1   ][Y(1) ]=[B*DX^(2)]
'THE SECCOND     [1  AA  1][Y(2) ]=[B*DX^(2)]
'THE LAST        [   1  AA][Y(19)]=[B*DX^(2)]
'THE EQ. SYSTEM IS SOLVED BY A GAUSS ROUTINE

     DIM C(80,80)
    '************* change M according to suggestion *******
     M=19   'Number of equations (9, 19, 39, 79)
    '*************     Study convergence            *******
'   INDATA
    A=.25
    B=8
    xin=0
    xend=10
    yin=0
    yend=0
    DX=(xend-xin)/(M+1)
    CC=B*DX^2
    AA=A*DX^(2)-2

'  SETUP OF COEFF. MATRIX
    n1=M
    for i=1 to n1
      C(i,n1+1)=CC
    next i

    for i=1 to n1-1
      C(i,i+1)=1
      C(i,i)=AA
    next i
      C(n1,n1-1)=1
      C(n1,n1)=AA

    for i=2 to n1-1
      C(i,i-1)=1
    next i
'   END OF SETUP

    print 0 ;"   "; yin
    gosub [gauss]
    print n1+1;"   ";yend
    print "analytical solution at x=5 ,10  20 and 40 71.9429"
   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)
     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
  NEXT J

     FOR J=1 TO M
   PRINT J;"   ";C(J,N)'solution vector
   REM
 NEXT J
return