Finite Differences 2D Not OK

MODULE Finite_differences
USE, intrinsic :: iso_fortran_env, only: real64
USE DATA
USE Aproximacion
USE Linear_System_Solver, ONLY: Band_Relaxation_Method
PRIVATE

PUBLIC:: Solver_BVP

CONTAINS

    SUBROUTINE Solver_BVP()
    !Ax=b
    REAL(kind=real64),DIMENSION(Extension, Extension)::Mat_Coeficientes !A
    REAL(kind=real64),DIMENSION(Extension)::Mat_Independiente !B
    REAL(kind=real64),DIMENSION(Extension):: Temperatura !x

    REAL(kind=real64):: Pos_x,Pos_y
    REAL(kind=real64):: Temperatura_Frontera
    REAL(kind=real64):: Fuente

    INTEGER I,J,Nodo_Iter

    OPEN(Unit=1,File="write_outs.dat")

    Mat_Coeficientes=0.0_real64; Mat_Independiente=0.0_real64; Temperatura=0.0_real64   !Inicializacion de las matrices en cero

    DO I=1,Extension                        !Asigno a toda la matriz de coeficientes la matriz identidad => Mat_Coeficientes=I
        Mat_Coeficientes(i,i)=1.0_real64
    END DO

!----------------------------ASIGNACION DE VALORES DE FRONTERA---------------------------------------------------------------

!-------------------------------FRONTERA DERECHA-----------------------------------------------------------------------------
    DO I= 1,Segmentos_y+1                           !Cantidad de nodos en el eje Vertical
        Pos_x = X_Final                             !Me mantengo siempre en el borde de la derecha
        Pos_y = Y_initial + (i-1)*Incremento_y      !Voy subiendo a medida que pasan las iteraciones
        CALL Condicion_Borde_Der(Pos_x,Pos_y,Temperatura_Frontera)
        Mat_Independiente(I*(Segmentos_x+1))=Temperatura_Frontera
    END DO

!-------------------------------FRONTERA IZQUIERDA----------------------------------------------------------------------------
    DO I= 1,Segmentos_y+1                               !Cantidad de nodos en el eje vertical
        Pos_x = X_Initial                               !Me mantengo siempre en el borde de la izquierda
        Pos_y = Y_initial + (i-1)*Incremento_y          !Voy subiendo a medida que pasan las iteraciones
        CALL Condicion_Borde_Izq(Pos_x,Pos_y,Temperatura_Frontera)
        Mat_Independiente((i-1)*(Segmentos_x+1)+1)=Temperatura_Frontera
    END DO

!------------------------------FRONTERA INFERIOR--------------------------------------------------------------------------------
    DO I= 1,Segmentos_x+1                               !Cantidad de nodos en el eje horizontal
        Pos_x = X_Initial + (i-1)*Incremento_x          !Voy desplazandome hacia la derecha a medida que pasan las iteraciones
        Pos_y = Y_initial                               !Me mantengo en la frontera inferior
        CALL Condicion_Borde_Inf(Pos_x,Pos_y,Temperatura_Frontera)
        Mat_Independiente(I)=Temperatura_Frontera
    END DO

!------------------------------FRONTERA SUPERIOR---------------------------------------------------------------------------------
    DO I= 1,Segmentos_x+1                               !Cantidad de nodos en el eje horizontal
        Pos_x = X_Initial + (i-1)*Incremento_x          !Voy desplazandome hacia la derecha a medida que pasan las iteraciones
        Pos_y = Y_Final                                 !Me mantengo en la frontera inferior
        CALL Condicion_Borde_Sup(Pos_x,Pos_y,Temperatura_Frontera)
        Mat_Independiente((Extension-Segmentos_x)+(I-1))=Temperatura_Frontera
    END DO

!--------------------------------------------------------------------------------------------------------------------------------

!----------------------------ASIGNACION DE VALORES DE MAT_INDEPENDIENTE (FUENTE)-------------------------------------------------

    DO I= 1,Segmentos_y                                 !Cantidad de nodos internos en el eje vertical

        DO J= 1,Segmentos_x                             !Cantidad de nodos internos en el eje horizontal

            Pos_x= X_Initial + j*Incremento_x           !Me voy desplazando horizontalmente hacia la derecha a medida que aumenta J
            Pos_y= Y_Initial + i*Incremento_y           !Me voy desplazando verticalmente hacia arriba a medida que aumenta I

            CALL Fuente_de_Calor(Pos_x,Pos_y,Fuente)

            Mat_Independiente(I*(Segmentos_x+1)+1+J) = Fuente/Conductividad

        END DO
    END DO

!----------------------------------------------------------------------------------------------------------------------------------

!-------------------------ASIGNACION DE VALORES MATRIZ DE COEFICIENTES-------------------------------------------------------------

    DO J = 1, Segmentos_y-1                 !Cantidad de veces que debo subir para cubrir todos los nodos incognita

        I = J*(Segmentos_x+1)               !Me ubica en el nodo final de cada fila

        DO Nodo_iter = I+2, I+Segmentos_x   !Recorro los nodos internos en cada fila

        Mat_Coeficientes(Nodo_Iter,Nodo_Iter) =-2.0_real64*(Cte_x+Cte_y)
        Mat_Coeficientes(Nodo_Iter,Nodo_Iter-1) = Cte_x
        Mat_Coeficientes(Nodo_Iter,Nodo_Iter+1) = Cte_x
        Mat_Coeficientes(Nodo_Iter,Nodo_Iter - Segmentos_x - 1) = Cte_y
        Mat_Coeficientes(Nodo_Iter,Nodo_Iter + Segmentos_x + 1) = Cte_y

        END DO

    END DO


    CALL Band_Relaxation_Method(Mat_Coeficientes,Mat_Independiente,Temperatura,Extension)

    DO I=1,Segmentos_y+1

        DO J=1,Segmentos_x+1

            Pos_x= X_Initial + (j-1)*Incremento_x           !Me voy desplazando horizontalmente hacia la derecha a medida que aumenta J
            Pos_y= Y_Initial + (i-1)*Incremento_y           !Me voy desplazando verticalmente hacia arriba a medida que aumenta I

        Write(1,*) Pos_x, Pos_y, Temperatura((I-1)*(Segmentos_x+1)+J)

        END DO

    END DO



END SUBROUTINE Solver_BVP

END MODULE Finite_Differences