module sistemaEc
    contains

    subroutine resolverSistema(sigmaA, muA)

        integer :: i, n, nmax, m
        real(8) :: errf, errx, eps
        real(8), dimension(2) :: x, y
        real(8), dimension(2,2) :: A
        real(8), intent(inout) :: sigmaA, muA

        !ITERACIONES MÁXIMAS, TOLERANCIA, Nº ECUACIONES
        nmax = 10
        eps = 1d-6
        m = 2

        !INICIALIZAMOS EL VECTOR SOLUCIÓN INICIAL, tal que x(:)=(/muA, sigmaA/)
        x = (2d0,6d0)

        !BUCLE PARA REALIZAR LAS APROXIMACIONES A LA SOLUCIÓN BUSCADA
        do n=1, nmax
            call f(x,y,m)
            call jacob_f(x,a,m)
            call factorizar (a,m)
            call sustituir (a,y,m)

            x = x + y

            errx = norma_2(y,m)/(max(eps,norma_2(x,m)))
            errf = norma_2(y,m)

            write(*,*) "ITERACION: ", n
            write(*,*) "sigmaA, muA: ", x(2), x(1)
            write(*,*)

            if (max(errx,errf) .lt. eps) then
                !GUARDAMOS LOS VALORES DE muA, sigmaA
                muA = x(1)
                sigmaA = x(2)
                write(*,*) "--------------"
                exit
            end if
        end do


        !SI NO CONVERGE, PARAMOS EL PROGRAMA
        if (n .gt. nmax) then
            write(*,'("numero maximo de iter.")')
            stop
        end if

        50 format (i2,4(1x,d13.7))

    end subroutine