'PROGRAM OPTIMER.BAS
 'PROGRAM WRITTEN BY INGEMAR BJERLE MARCH 2002
 'FINDING OPTIMAL TEMPERATURE PROFILE IN BATCH REACTOR
 'REACTIONS:   A ---> B ---> C
 '                K1     K2
 '        DCA/DT = -K1*A^2   DCB/DT= K1*A^2 - K2*B
 'MAXIMATION OF THE AMOUNT OF B
 'THE EQUATIONS DEDUCED BY VARIATIONAL CALCULUS

    INPUT "tmax (TRY 400)=";T1
    INPUT "tmin (TRY 300)=";T2
    INPUT "tstart(300-400=";T3
    PRINT""
    print "         ITERATIONS"
    print"        A             B"
    F1=4000:F2=620000:F3=.503271
    DIM A(101):DIM B(101):DIM T(101):DIM K1(101):DIM K2(101)
    FOR I=0 TO 100: T(I)=T3: NEXT I
    QQ=0

    [CALC]
    T=0: D=.01: A1=1: A(0)=1: B1=0: T(0)=T1: K3=-5000/1.987: E=5*10^7
    FOR I=0 TO 100   'EULER FORWARD
      K1(I)=F1*EXP(K3/T(I))
      K2(I)=F2*K1(I)*K1(I)/F1/F1
      D1= (-1)* K1(I)*A1*A1*D
      B1= B1-D1-K2(I)*B1*D: A1=D1+A1: T=T+D
      A(I+1)=A1: B(I+1)=B1
    NEXT I

     QQ=QQ+1
     PRINT QQ;"  "; A(100);"        ";B(100)
     IF ABS(( B(100)-Q)/ B(100))<.0000001 GOTO [PRINT]
     Q=B(100)

     L1=0: L2=1: T=1
     FOR I=100 TO 0 STEP -1'  EULER BACKWARD
       S=(L2-L1)*K1(I)*A(I)*A(I)/2: M=L2*K2(I)
       T(I)=T(I)+E*F3/T(I)/T(I)*(S-M*B(I))
       L1=L1+S/A(I)*4*D
       L2=L2-M*D
       IF T(I)>T1 THEN T(I)=T1
       IF T(I)<T2 THEN T(I)=T2
     NEXT I

    GOTO [CALC]

    [PRINT] print " "
     print"        optimal temperature profile"
     print"         time        temperature"
     T=0
     FOR I=0 TO 10
       PRINT USING("##########.###",T);USING("############.##",T(I))
       T=T+1*D
     NEXT I
     T=.2
     PRINT USING("##########.###",T);USING("############.##",T(20))
     FOR I=40 TO 100 STEP 20
       T=T+20*D
       PRINT USING("##########.###",T);USING("############.##",T(I))
     NEXT I
     PRINT "PRESS ENTER TO END"
     INPUT A$ 
     PRINT"                    *** end ***"
     END