How to compute the Jacobian of a Vector-to-Vector complicated function

I am working on a modified finite element simulation. I have a complicated function to project a 648 units long vector to a 648 units long vector (i.e., from nodal displacements to nodal forces, the magnitude will be maginified to 10000-to-10000). The function (project) is kind of complicated, I put it as follows. My goal is to get the Jacobian matirix of the Vector-to-Vector function at every Vector input. I have tried some automatic differentiation packages, as well as Zygote, but failed always. Every effective recommendation or clue is appreciated.

I also upload the data (the constants in the fucntion) needed as a jld2 in my Github GitHub - DonggeJia/the-function-needs-Jacobian: A project from nodal displacement to nodal force of LDPM

using PhysicalMeshes
function constitutive(noddis)

    #get volumetric strain for every element
    V_epsi = vec(zeros(length1, 1))
    division = vec(zeros(length1, 1))

    newverticecoor = verticeoftetrahedron1 + [noddis[(jj*6-5):(jj*6-3)] for jj = 1:length_nod][tttt1]
    for tnumber = 1:length2
        epsiVforte = (PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(newverticecoor[tnumber, 1]), PVector(newverticecoor[tnumber, 2]), PVector(newverticecoor[tnumber, 3]), PVector(newverticecoor[tnumber, 4]))) -
                      PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))) / PhysicalMeshes.volume(PhysicalMeshes.Tetrahedron(PVector(verticeoftetrahedron1[tnumber, 1]), PVector(verticeoftetrahedron1[tnumber, 2]), PVector(verticeoftetrahedron1[tnumber, 3]), PVector(verticeoftetrahedron1[tnumber, 4])))

        for numele = 1:length1
            if elements41[numele][1] ∈ tetrahedron1[tnumber, :] && elements41[numele][2] ∈ tetrahedron1[tnumber, :]
                V_epsi = V_epsi + vcat(vec(zeros(numele-1,1)), epsiVforte, vec(zeros(length1 - numele,1)))
                division = division + vcat(vec(zeros(numele-1,1)), 1, vec(zeros(length1 - numele,1)))
                

            end
        end
    end
    V_epsi_1 = deepcopy(@. V_epsi / division)
    #get other strains for every element

    nodalforce1 = vec(zeros(length(pointsfinal)*6,1))
    for kjh = 1:length1

        epsiNk = BNjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BNik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiMk = BMjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BMik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiLk = BLjk2[kjh] * noddis[(6*elements41[kjh][2]-5):(6*elements41[kjh][2])] - BLik2[kjh] * noddis[(6*elements41[kjh][1]-5):(6*elements41[kjh][1])]
        epsiTk = sqrt(epsiMk^2 + epsiLk^2)

        if epsiNk >= 0
            effectiveepsi = sqrt(epsiTk^2 * alpha + epsiNk^2)
            # tensile-shear boundary
            omig = atan(epsiNk / sqrt(alpha) / epsiTk)
            s = sin(omig)
            c = cos(omig)
            if omig == atan(sqrt(alpha / mu))
                effecstress0b = 0.5 * (sigmt + 2 * sigma) * sigmt / ((sigmt + sigma) * sin(omig))
            else
                effecstress0b = (-(sigmt + sigma) * s + sqrt(((sigmt + sigma) * s)^2 + (alpha * (c / mu)^2 - s^2) * (sigmt + 2 * sigma) * sigmt)) / (alpha * (c / mu)^2 - s^2)
            end
            Hs = 2 * alpha * EN / ((2 * alpha * EN * Gs) / (sigmt^2 * le2[kjh]) - 1)
            Ht = 2 * EN / ((2 * EN * Gt) / (sigmt^2 * le2[kjh]) - 1)
            H = Hs + (Ht - Hs) * (2 * omig / pi)^nt

            effectstressb = effecstress0b * exp(-H / effecstress0b * maximum([maxeffecsigm[kjh] - effecstress0b / EN, 0]))

            if effectiveepsi * EN <= effectstressb
                effectivestress = effectiveepsi * EN
            else
                effectivestress = effectstressb
            end
            if epsiNk == 0.0
                sigmNk = 0.0
            else
                sigmNk = effectivestress / effectiveepsi * epsiNk
            end

            sigmTk = effectivestress / effectiveepsi * epsiTk * alpha

            if epsiLk == 0.0
                sigmLk = 0.0
            else
                sigmLk = sigmTk * epsiLk / epsiTk  ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk
            end
            if epsiMk == 0.0
                sigmMk = 0.0
            else
                sigmMk = sigmTk * epsiMk / epsiTk
            end

        else #epsiNk < 0
            #compressive
            rdv = (epsiNk - V_epsi_1[kjh]) / V_epsi_1[kjh]
            if V_epsi_1[kjh] >= 0
                if epsiNk >= sigmc0 / EN
                    sigmNk = epsiNk * EN
                else
                    sigmNk = sigmc0
                end
            else #V_epsi_1[kjh] < 0
                H = Hc / (1 + kapha * rdv)
                sigmcb = sigmc0 * exp(H / sigmc0 * maximum([epsic0 - V_epsi_1[kjh], 0]))
                if epsiNk >= sigmcb / EN
                    sigmNk = epsiNk * EN
                else
                    sigmNk = sigmcb
                end
            end
            #shear stress under compression

            if sigmNk <= sigmN0
                omig = atan(epsiNk / sqrt(alpha) / epsiTk)
                s = sin(omig)
                c = cos(omig)
                sigmTb0 = -sigmc0 / sqrt(s^2 + alpha * c^2 / betta)
            else
                sigmTb0 = mu * sqrt((epsiNk - sigmt - sigma)^2 - sigma^2)
            end
            epsiT0 = sigmTb0 / ET

            omig = atan(epsiNk / sqrt(alpha) / epsiTk)
            Hs = 2 * alpha * EN / ((2 * alpha * EN * Gs) / (sigmt^2 * le2[kjh]) - 1)
            H = Hs - Hs * (-2 * omig / pi)^nc
            sigmTb = sigmTb0 * exp(-H / sigmTb0 * maximum([epsiTk - epsiT0, 0]))

            if epsiTk < sigmTb / ET
                sigmTk = epsiTk * ET
            else
                sigmTk = sigmTb
            end
            if epsiLk == 0.0
                sigmLk = 0.0
            else
                sigmLk = sigmTk * epsiLk / epsiTk ##the ratio between L and M, use epsiLk / epsiTk maintains the relation between lk and mk
            end
            if epsiMk == 0.0
                sigmMk = 0.0
            else
                sigmMk = sigmTk * epsiMk / epsiTk
            end

        end

        nodalforce1 = nodalforce1 + vcat(vec(zeros(6*elements41[kjh][1]-6,1)), (-le2[kjh] * Ak2[kjh] * (sigmNk * BNik2[kjh] + sigmMk * BMik2[kjh] + sigmLk * BLik2[kjh]))',vec(zeros(length_nod*6 - (6*elements41[kjh][1]),1)))
        nodalforce1 = nodalforce1 +[vec(zeros(1,6*elements41[kjh][2]-6))' (-le2[kjh] * Ak2[kjh] * (sigmNk * BNjk2[kjh] + sigmMk * BMjk2[kjh] + sigmLk * BLjk2[kjh])) vec(zeros(1,length_nod*6 - (6*elements41[kjh][2])))']'
    end
    return nodalforce1
end

You could try a brute force finite difference Jacobian. If your Jacobian is sparse, look into SparseDiffTools.jl otherwise it’s not hard to code this yourself (see Chapters 1 and 2 of this suite of notebooks).

Thank you so much, I am trying.