Off[General::"spell1"];
PREP3 := (eg2[u_] := -4*
              NIntegrate[
                BesselJ[1, z] BesselJ[0, z]*
                  Exp[-u*2*z]/(z*(1 + Exp[-u*2*z])), {z, 0, 
                  Infinity}]; 
F[x_] := 
            Sin[x]*Sqrt[1 + (Ud/Cos[x])^2]; 
Fd[
              x_] := (Ud^2 + Cos[x]^4)/(Cos[x]^2*
                  Sqrt[Ud^2 + (Cos[x])^2]); 
xl = 
            Table[((j -  .5)/2/LL)*Pi, {j, 1, LL}]; tx = F[xl]; 
          pc = Pi/2/Ud; 
Print[" tx=", 
            tx]; 
gg[x_, y_] := 
            Cosh[pc*x]*Fd[y]* .5*
              Sech[pc*F[
                      y]]/(1 + (Sinh[pc*x]*
                          Sech[pc*F[y]])^2); 
ss = 
            Table[Table[
                NIntegrate[
                    gg[tx[[i]], y], {y, (j - 1)*Pi/(2*LL), 
                      j*Pi/(2*LL)*(1.  -  .00001*KroneckerDelta[i, LL]*
                              KroneckerDelta[j, LL])}, 
                    MaxRecursion -> 20]/Ud, {j, LL}], {i, 
                LL}]; 
Print[
            Table[Sum[ss[[i]][[j]], {j, LL}], {i, 
                LL}]]; 
kk = 
            Table[Sin[Pi*(i -  .5)/(2*LL)], {i, 
                LL}]; 
s0 = 
            Table[NIntegrate[
                NIntegrate[ .5/Ud*
                      Sech[Pi (F[y] - Sin[z1])/(2*
                                  Ud)]/(2 Pi), {z1, 0, 2*Pi}]*
                  Fd[y], {y, (j - 1)*Pi/(2*LL), 
                  j*Pi/(2*LL)}], {j, 
                LL}]; 
Print[" s0=", s0]; 
          kx = Table[
              Pi*(i -  .5)/LL, {i, 
                LL}]; 
Print[" kx= ", 
            kx]; 
rou = 
            Table[(1.  + 
                    Cos[kx[[i]]]*
                      NIntegrate[ .5*
                          Re[1. /
                              Sqrt[1 - (Sin[kx[[i]]] - z*Ud - 
                                        I*Ud)^2]]*
                          Sech[Pi*z/2], {z, -Infinity, Infinity}])/
                LL, {i, LL}]; 
Print[" ro*2*Pi= ", 
            rou*LL]; 
sk = 
            Table[Table[
                NIntegrate[
                    gg[Sin[kx[[i]]], y], {y, (j - 1)*Pi/(2*LL), 
                      j*Pi/(2*LL)}, MaxRecursion -> 20]/Ud, {j, 
                  LL}], {i, LL}]; 
slk = 
            Table[Table[(Sech[(tx[[i]] - Sin[kx[[j]]])*
                          Pi/(2*Ud)] + 
                      Sech[(tx[[i]] + Sin[kx[[j]]])*
                          Pi/(2*Ud)])*Cos[kx[[j]]]*
                  Pi/(4*LL*Ud), {j, LL}], {i, 
                LL}]; 
re = 
            Re[Sqrt[1 - (tx - I*Ud)^2]]; 
          k0 = -2*Cos[kx] - 
              4*sk . re; 
Print[" k0= ", k0]; 
          e0 = -2 Cos[kx] . rou; 
Print[" e0= ", 
            e0, "     e-exact=", eg2[Ud]]; 
s1 = 
            Table[slk[[i]] . k0, {i, 
                LL}]; 
Print[" s1= ", 
            s1]);
\(\(XXX6 := \((LC21 = 2*LC + 1; 
dx1 = 
            Table[Table[0. , {j, LL}], {i, LC21}]; 
          AO2 = Table[\(-ss\), {j, LC*2}]; 
AO2[\([LC]\)] = 
            slk; AO2[\([LC + 1]\)] = sk; 
          AO3 = Table[\(-ss\), {j, 
                LC*2}]; 
AO3[\([LC]\)] = \(-sk\); 
          AO3[\([LC + 1]\)] = 
            slk; 
Linear2 := \((AA1 = 
                Table[Table[
                    Table[KroneckerDelta[i, j], {i, LL}], {j, LL}], {k, 
                    LC21}]; 
AA1[\([1]\)] = 
                AA1[\([1]\)] - 
                  Table[ss[\([j]\)]*ffb2, {j, 
                      LL}]; 
AA1[\([LC21]\)] = 
                AA1[\([LC21]\)] - 
                  Table[ss[\([j]\)]*ffa2, {j, LL}]; 
AA2 = 
                Table[Table[\(AO2[\([j]\)]\)[\([k]\)]*\((1. /
                          df[\([j + 1]\)])\), {k, LL}], {j, 
                    LC*2}]; 
AA3 = 
                Table[Table[\(AO3[\([j]\)]\)[\([k]\)]*\((1. /
                          df[\([j]\)])\), {k, LL}], {j, 
                    LC*2}]; 
WW = Inverse[AA1[\([1]\)]]; 
              AA2[\([1]\)] = 
                WW . AA2[\([1]\)]; 
dinx[\([1]\)] = 
                WW . dinx[\([1]\)]; 
Do[{AA1[\([j]\)] = 
                    AA1[\([j]\)] - 
                      AA3[\([j - 1]\)] . 
                        AA2[\([j - 1]\)]; 
dinx[\([j]\)] = 
                    dinx[\([j]\)] - 
                      AA3[\([j - 1]\)] . 
                        dinx[\([j - 1]\)]; 
WW = 
                    Inverse[AA1[\([j]\)]]; 
                  dinx[\([j]\)] = WW . dinx[\([j]\)]; 
If[
                    j < LC21, AA2[\([j]\)] = WW . AA2[\([j]\)]]}, {j, 2, 
                  LC21}]; 
Do[{dinx[\([j]\)] = 
                    dinx[\([j]\)] - AA2[\([j]\)] . dinx[\([j + 1]\)]}, {j, 
                  2*LC, 1, \(-1\)}])\); 
If[
            inicnt == 1, {dl = dx1; 
              dl[\([LC + 1]\)] = 
                Table[Log[Cosh[mb]/Cosh[hb]], {i, 
                    LL}]; 
Do[{dl[\([LC + 1 + i]\)] = 
                    Table[Log[\((Sinh[\((i + 1)\) \((hb +  .00001)\)]/
                                Sinh[hb +  .00001])\)^2 - 1. ], {j, 
                        LL}]; 
dl[\([LC + 1 - i]\)] = 
                    Table[Log[\((Sinh[\((i + 1)\) \((mb +  .00001)\)]/
                                Sinh[mb +  .00001])\)^2 - 1. ], {j, 
                        LL}]}, {i, LC}]}]; 
dl2 = dx1; 
          dlp1 = Log[Exp[dl] + 1]; 
q = 0; 
          prr = {" Newton Method"}; 
Label[
            beg]; 
dl2[\([LC + 1]\)] = 
            k0*bet + 
              sk . \((dlp1[\([LC]\)] - 
                    dlp1[\([LC + 2]\)])\); 
dl2[\([LC + 
                  2]\)] = 
            bet*s1 - slk . dlp1[\([LC + 1]\)] + 
              ss . dlp1[\([LC + 
                      3]\)]; 
dl2[\([LC]\)] = \(-slk . 
                  dlp1[\([LC + 1]\)]\) + 
              ss . dlp1[\([LC - 1]\)]; 
dl2[\([LC21]\)] = 
            ss . \((dlp1[\([2*LC]\)] + 
                  2. *Log[
                      Exp[ .5*dlp1[\([LC21]\)]]*Cosh[hb] + 
                        Sqrt[1.  + \((Sinh[hb]*
                                  Exp[ .5*
                                      dlp1[\([LC21]\)]])\)^2]])\); 
dl2[\([1]\)] = 
            ss . \((dlp1[\([2]\)] + 
                  2. *Log[
                      Exp[ .5*dlp1[\([1]\)]]*Cosh[mb] + 
                        Sqrt[1.  + \((Sinh[mb]*
                                  Exp[dlp1[\([1]\)]* .5])\)^2]])\); 
Do[{dl2[\([LC + 1 + j]\)] = 
                ss . \((dlp1[\([LC + j]\)] + 
                      dlp1[\([LC + 2 + 
                            j]\)])\); 
dl2[\([LC + 1 - 
                      j]\)] = 
                ss . \((dlp1[\([LC - j]\)] + dlp1[\([LC + 2 - j]\)])\)}, {j, 
              2, LC - 1}]; 
dinx = 
            dl2 - dl; 
err2 = 
            Sum[Sum[\(\(Abs[dinx]\)[\([i]\)]\)[\([j]\)], {j, LL}], {i, 
                LC21}]; 
df = 
            1.  + Exp[\(-dl\)]; 
ffa2 = \((1.  - 
                Cosh[hb]/
                  Sqrt[\((1 + Exp[dl[\([LC21]\)]])\) \((1 + 
                          Sinh[hb]^2 \((1 + 
                                Exp[dl[\([LC21]\)]])\))\)])\); 
ffb2 = \((1.  - 
                Cosh[mb]/
                  Sqrt[\((1 + Exp[dl[\([1]\)]])\) \((1 + 
                          Sinh[mb]^2 \((1.  + 
                                Exp[dl[\([1]\)]])\))\)])\); 
Linear2; 
dl = dl + dinx; 
          dlp1 = Log[1.  + Exp[dl]]; 
q += 
            1; 
prr = Append[prr, {q, err2}]; 
          If[q < 10 && err2 > 10. ^\((\(-11\))\), 
            Goto[beg]]; 
gp = 
            e0 - \((rou . dlp1[\([LC + 1]\)] + s0 . dlp1[\([LC + 2]\)])\)/
                bet - 2*Ud + 
              mb/bet; 
 (*\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(----\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)*) din1 = dx1; 
          din1[\([LC + 1]\)] = k0; 
          din1[\([LC + 2]\)] = s1; 
din2 = 
            dx1; 
din2[\([1]\)] = 
            2*ss . \((Sinh[mb]/
                    Sqrt[Sinh[mb]^2 + 
                        Exp[\(-dlp1[\([1]\)]\)]])\); 
din3 = dx1; 
din3[\([LC21]\)] = 
            2*ss . \((Sinh[hb]/
                    Sqrt[Sinh[hb]^2 + 
                        Exp[\(-dlp1[\([LC21]\)]\)]])\); 
dinx = 
            Table[Table[{\(din1[\([i]\)]\)[\([j]\)], \(din2[\([i]\)]\)[\([j]\)], \(din3[\([i]\)]\)[\([j]\)]}, {j, LL}], {i, 
                LC21}]; 
          Linear2; 
dd = 
            rou . \((dinx[\([LC + 1]\)]/df[\([LC + 1]\)])\) + 
              s0 . \((dinx[\([LC + 2]\)]/
                    df[\([LC + 2]\)])\); 
eng2 = 
            e0 - dd[\([1]\)] - 2*Ud*dd[\([2]\)]; eld2 = 1.  - dd[\([2]\)]; 
          mag2 = dd[\([3]\)]; 
 (*\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(--\(----\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)\)*) ent = 
            bet*\((eng2 - gp)\) - 
              hb*mag2 + \((mb - 2*Ud*bet)\)*eld2; 
diny = 
            dinx; 
dinx = 
            Table[Table[{\(\(diny[\([i]\)]\)[\([j]\)]\)[\([1]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([1]\)], \(\(diny[\([i]\)]\)[\([j]\)]\)[\([1]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([2]\)], \(\(diny[\([i]\)]\)[\([j]\)]\)[\([1]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([3]\)], \(\(diny[\([i]\)]\)[\([j]\)]\)[\([2]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([2]\)], \(\(diny[\([i]\)]\)[\([j]\)]\)[\([2]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([3]\)], \(\(diny[\([i]\)]\)[\([j]\)]\)[\([3]\)]*\(\(diny[\([i]\)]\)[\([j]\)]\)[\([3]\)]}, {j, LL}], {i, LC21}]; 
dy1 = dinx[\([1]\)]; 
          dy2 = dinx[\([LC21]\)]; 
          dinx = dinx*\((1.  - 1. /df)\); 
dz = 
            dlp1[\([LC21]\)]; z1 = Exp[dz]; 
          z2 = 1 + Sinh[hb]^2*z1; 
dz = dlp1[\([1]\)]; 
          z1 = Exp[dz]; 
          z2 = 1 + Sinh[mb]^2*z1; 
g1 = 
            2*Cosh[mb]*\((Exp[ .5*
                      dlp1[\([1]\)]]/\((Sqrt[
                        z2])\)^3)\); 
g2 = 
            Sinh[mb]*
              Exp[dl[\([1]\)]]/Sqrt[z1*z2^3]; 
g3 = \(-1\)/
                z1 + Cosh[mb]/
                  Sqrt[z1*z2]*\((1 - 1/z2 + 
                    1.5/\((z1*z2)\))\); 
dz = 
            dlp1[\([LC21]\)]; z1 = Exp[dz]; 
          z2 = 1 + Sinh[hb]^2*z1; 
f1 = 
            2*Cosh[hb]*\((Exp[ .5*
                      dlp1[\([LC21]\)]]/\((Sqrt[
                        z2])\)^3)\); 
f2 = 
            Sinh[hb]*
              Exp[dl[\([LC21]\)]]/
                Sqrt[z1*z2^3]; 
f3 = \(-1\)/z1 + 
              Cosh[hb]/
                  Sqrt[z1*z2]*\((1 - 1/z2 + 
                    1.5/\((z1*z2)\))\); 
dy1 = dy1*g3; 
          dy2 = dy2*f3; 
          ll1 = Transpose[diny[\([1]\)]]; 
ll2 = 
            Transpose[diny[\([LC21]\)]]; 
dz1 = 
            Transpose[dy1]; 
          dz1[\([4]\)] = 
            dz1[\([4]\)] + g1 + 
              2*g2*ll1[\([2]\)]; 
dz1[\([2]\)] = 
            dz1[\([2]\)] + g2*ll1[\([1]\)]; 
          dz1[\([5]\)] = 
            dz1[\([5]\)] + g2*ll1[\([1]\)]; 
dy1 = 
            ss . Transpose[dz1]; 
dz1 = Transpose[dy2]; 
          dz1[\([6]\)] = 
            dz1[\([6]\)] + f1 + 
              2*f2*ll2[\([3]\)]; 
dz1[\([2]\)] = 
            dz1[\([2]\)] + f2*ll2[\([1]\)]; 
          dz1[\([5]\)] = 
            dz1[\([5]\)] + f2*ll2[\([1]\)]; 
dy2 = 
            ss . Transpose[dz1]; 
dinx[\([1]\)] = 
            dinx[\([1]\)] + dy1; dinx[\([LC21]\)] = dinx[\([LC21]\)] + dy2; 
          Linear2; 
dd2 = 
            rou . \((dinx[\([LC + 1]\)]/df[\([LC + 1]\)])\) + 
              s0 . \((dinx[\([LC + 2]\)]/
                    df[\([LC + 
                          2]\)])\); 
u1 = \((dd2[\([3]\)]*
                    dd2[\([4]\)] - 
                  dd2[\([2]\)]*dd2[\([5]\)])\)/\((dd2[\([5]\)]^2 - 
                  dd2[\([4]\)]*
                    dd2[\([6]\)])\); 
u2 = \((dd2[\([2]\)]*
                    dd2[\([6]\)] - 
                  dd2[\([3]\)]*dd2[\([5]\)])\)/\((dd2[\([5]\)]^2 - 
                  dd2[\([4]\)]*
                    dd2[\([6]\)])\); 
 (*\(--\(-specific\)\)\ heat\ in\ constant\ density\ and\ magnetization, 
            susceptibility\ in\ constant\ density\ and\ temperature, \(\(compressibility\ in\ constant\ temperature\ and\ mag/\(eld--\)\)\(-\)\)*) 
            sph = \((dd2[\([1]\)] + u1*dd2[\([3]\)] + u2*dd2[\([2]\)])\)*
              bet^2; 
chi = 
            bet*\((dd2[\([6]\)] - 
                  dd2[\([5]\)]^2/dd2[\([4]\)])\); 
comp = 
            bet*\((dd2[\([4]\)] - \((dd2[\([3]\)]*eld2 + dd2[\([4]\)]*mag2)\)*
                    dd2[\([5]\)]/\((dd2[\([6]\)]*eld2 + 
                          dd2[\([3]\)]*
                            mag2)\))\); 
PRI := \((Print["U=  ", 4*Ud, ", LL=", LL, ", LC=", LC, ", Temp=", 
                1. /bet, ", beta=", bet, ",  U/2-A=", 
                mb/bet, ",  h=", hb/bet, "---------------------"]; 
              Print[prr]; 
              Print[dd, 
                dd2]; 
Print[" Grand Potential=", 
                InputForm[gp], ", Internal Energy=", 
                InputForm[eng2], ", Electron Density=", 
                InputForm[eld2], ", Magnetization=", 
                InputForm[mag2], ", Entropy=", 
                InputForm[
                  ent]]; 
Print[" specific heat=", 
                InputForm[sph], ", susceptibility=", 
                InputForm[chi], ", compressibility=", 
                InputForm[comp]])\); 
If[prcont == 1, 
            PRI])\);\)\)
Search := (q2 = 0; prcont = 0; Label[beg2]; XXX6; 
        q2 += 1; 
Print[q2, ",  ", q, ",  ", 
          mb/bet, ",     ", eld2 - tgd, ",     ", 
          eld2]; 
mb = 
          Min[Max[mb + (eld2 - tgd)/Abs[dd2[[4]]], mb/2], 2*mb]; 
        inicnt = 0; 
If[
          Abs[eld2 - tgd] > 10. ^(-8) && q2 < 10, Goto[beg2]]; 
        PRI)

$DefaultFont = {"Bold", 20};
Print[Timing[LL = 16; Ud = 1. ; PREP3]]; LC = 6; mb = 0; hb = 0; 
    inicnt = 1; prcont = 1; Do[{bet =  .01*j; XXX6}, {j, 10}]

Print[Timing[LL = 16; Ud = 2. ; PREP3]]; LC = 6; mb = 0; hb = 0; 
    inicnt = 1; prcont = 1; Do[{bet =  .01*j; XXX6}, {j, 10}]