REM Resolucion numerica mediante  metodo EXPLICITO
REM de la EDP parabolica:
REM             du/dt=k d^2u/dx^2
REM con  condiciones de contorno u(0,t)=u(1,t)=0 y
REM condicion inicial u(x,0)=fnexacta(x,0)
REM        VER PAG. 162-163 DE KOONIN


CLS
PRINT "        METODO EXPLICITO"
PRINT
n = 30  'numero de puntos (escojase par)
ix = 1 / n
dt = .001      ' dt=constante difusiva*discretizado temporal
itmax = 50     ' numero maximo de iteraciones
saltot = 10    ' mostraremos los resultados cada saltot iteraciones
s = dt / ix ^ 2
PRINT "numero de puntos="; n
PRINT "inc_x="; ix;
PRINT "      k*inc_t="; dt
PRINT USING "dt/ix^2=##.#### <1/2 ==> estabilidad"; s
PRINT
DIM u(n)
DEF fngauss (x, t) = EXP(-20 * (x - .5) ^ 2 / (1 + 80 * t)) / SQR(1 + 80 * t)
DEF fnexacta (x, t) = fngauss(x, t) - fngauss(x - 1, t) - fngauss(x + 1, t)

t = 0
u(0) = 0: u(n) = 0
FOR i = 1 TO n - 1
        u(i) = fnexacta(i * ix, 0)
NEXT i
PRINT "Errores en el punto medio para varios tiempos:"
FOR it = 1 TO itmax
        uvieja = 0
        FOR i = 1 TO n - 1
            unueva = u(i) + s * (uvieja + u(i + 1) - 2 * u(i))
            uvieja = u(i)
            u(i) = unueva
        NEXT i
        IF it MOD saltot <> 0 THEN GOTO 100
        t = dt * it
        dif = fnexacta(ix * n / 2, t) - u(n / 2)
       PRINT USING "tiempo=#.######, error=#.#########"; t; dif
        
100 NEXT it

        PRINT
        PRINT USING "Valores en varias posiciones en t=#.######"; t
        FOR i = 1 TO n - 1
            IF i MOD 5 <> 0 THEN GOTO 200
            exacto = fnexacta(i * ix, t)
        PRINT USING "u(##)=#.####^^^^  exacta=#.####^^^^"; i; u(i); exacto
200     NEXT i