I think there is something wrong in this Mathematica Inverse Kinematic algorithm, because if change the sign of the fith joint alpha denavit hartenberg parameter the result is correct. The code is the following:
InverseKinematics6DOF[pe_, quat_, params_] := Module[
   {
    (*Extract position components*) 
    px, py, pz,
    (*Extract quaternion components*) 
    (*q0,q1,q2,q3,*)
    qw, qx, qy, qz,
    (*Compute rotation matrix from quaternion*) 
    R,
    (*End-effector orientation vectors*) 
    ne, se, ae,
    (*Wrist position*) 
    pw,
    (*Joint angles*) 
    q1Sol, q2Sol, q3Sol, q4Sol, q5Sol, q6Sol,
    (*Intermediate calculations*) 
    L2, L3, L4, A2, Pw,
    R03,
    T01,
    phi, beta, gamma, eta, r, a1 = a1 /. params, a2 = a2 /. params, 
    a3 = a3 /. params, d1 = d1 /. params, d4 = d4 /. params, 
    d6 = d6 /. params,
    (*List of joint angle solutions*) 
    jointAngles123, q456Sol, allSolutions
    }, {px, py, pz} = pe;
   {qw, qx, qy, qz} = quat;
   R03 = 
    Evaluate[(A0iall[[4]][[1 ;; 3, 1 ;; 3]]) /. params /. 
        q1 -> q1temp /. q2 -> q2temp /. q3 -> q3temp];
   T01 = Evaluate[A0iall[[2]] /. params /. q1 -> q1temp];
   (*Define your robot's DH parameters here*) 
   (* a1=a1/. params;a2=a2/. params;a3=a3/. params;
   d1=d1/. params;d4=d4/. params;d6=d6/. params;*)
   (*Step 2:Compute the Rotation Matrix from the Quaternion*) 
   R = {{1 - 2  qy^2 - 2  qz^2, 2  qx  qy - 2  qz  qw, 
      2  qx  qz + 2  qy  qw}, {2  qx  qy + 2  qz  qw, 
      1 - 2  qx^2 - 2  qz^2, 
      2  qy  qz - 2  qx  qw}, {2  qx  qz - 2  qy  qw, 
      2  qy  qz + 2  qx  qw, 1 - 2  qx^2 - 2  qy^2}};
   (*R={{1-2 (q2^2+q3^2),2 (q1 q2-q0 q3),2 (q1 q3+q0 q2)},{2 (q1 q2+
   q0 q3),1-2 (q1^2+q3^2),2 (q2 q3-q0 q1)},{2 (q1 q3-q0 q2),2 (q2 q3+
   q0 q1),1-2 (q1^2+q2^2)}};*)
   (*Step 3:Compute the End-Effector Orientation Vectors*)
    ne = R[[All, 1]];
   se = R[[All, 2]];
   ae = R[[All, 3]];
   (*Step 4:Compute the Wrist Position*) 
   pw = {px, py, pz} - d6*ae;
   Pw = Append[pw, 1];
   (*Step 5:Compute q1*) 
   (*q1Sol=ArcTan[pw[[2]],pw[[1]]];*)
   q1Sol = ArcTan[pw[[1]], pw[[2]]];
   q1Sol = {q1Sol, q1Sol + Pi};
   (*Two possible solutions*) 
   (*Step 6:Compute q2 and q3 for each q1*)
    L2 = a2;
   L3 = Sqrt[d4^2 + a3^2];
   A2 = a3;
   L4 = Sqrt[A2^2 + d4^2];
   phi = ArcCos[(A2^2 + L4^2 - L3^2)/(2 A2 L4)];
   (*Initialize list of joint angles*) 
   jointAngles123 = {};
   Do[
    (*Transform pw to frame 1*) 
    T01 = T01 /. q1temp -> q1p;
    pw1 = Inverse[T01] . Pw;
    (*pw1=RotationMatrix[-q1p,{0,0,1}].(pw-{a1,0,d1});*)
    r = Sqrt[pw1[[1]]^2 + pw1[[2]]^2];
    (*beta=ArcTan[-pw1[[2]],pw1[[1]]];*)
    beta = ArcTan[pw1[[1]], -pw1[[2]]];
    gamma = ArcCos[(L2^2 + r^2 - L3^2)/(2 L2 r)];
    q2ElbowUp = (Pi/2) - beta - gamma;
    q2ElbowDown = (Pi/2) - beta + gamma;
    eta = ArcCos[(L2^2 + L4^2 - r^2)/(2 L2 L4)];
    q3ElbowUp = Pi - phi - eta;
    q3ElbowDown = Pi - phi + eta;
    (*Append both elbow up and elbow down solutions*)
     AppendTo[jointAngles123, {q1p, q2ElbowUp, q3ElbowUp}];
    AppendTo[jointAngles123, {q1p, q2ElbowDown, q3ElbowDown}];, {q1p, 
     q1Sol}];
   (*Step 7:Compute q4,q5,q6 for each set of q1,q2,q3*) 
   q456Sol = 
    Table[Module[{R01, R12, R23, R36, q4Up, q5Up, q6Up, q4Down, 
       q5Down, q6Down, a33}, {q1sol, q2sol, q3sol} = joint;
      (*Compute the rotation matrix R0_3*)
      (* R01=RotationMatrix[q1,{0,0,1}];
      R12=RotationMatrix[q2,{0,1,0}];
      R23=RotationMatrix[q3,{0,1,0}];
      R03=R01.R12.R23;*)
      R03 = 
       R03 /. q1temp -> q1sol /. q2temp -> q2sol /. q3temp -> q3sol;
      R36 = Transpose[R03] . R;
      a33 = R36[[3, 3]];
      (*ifPositiveQ5=ArcTan[Sqrt[R36[[1,3]]^2+R36[[2,3]]^2],a33];
      ifNegativeQ5=ArcTan[-Sqrt[R36[[1,3]]^2+R36[[2,3]]^2],a33];*)
      ifPositiveQ5 = 
       ArcTan[a33, Sqrt[R36[[1, 3]]^2 + R36[[2, 3]]^2]];
      ifNegativeQ5 = ArcTan[a33, -Sqrt[R36[[1, 3]]^2 + R36[[2, 3]]^2]];
      (*Wrist up solution*)
      q5Up = ifPositiveQ5;
      (*q4Up=ArcTan[R36[[2,3]],R36[[1,3]]];
      q6Up=ArcTan[R36[[3,2]],-R36[[3,1]]];*)
      q4Up = ArcTan[R36[[1, 3]], R36[[2, 3]]];
      q6Up = ArcTan[-R36[[3, 1]], R36[[3, 2]]];
      (*Wrist down solution*) 
      q5Down = ifNegativeQ5;
      (*q4Down=ArcTan[-R36[[2,3]],-R36[[1,3]]];
      q6Down=ArcTan[-R36[[3,2]],R36[[3,1]]];*)
      q4Down = ArcTan[-R36[[1, 3]], -R36[[2, 3]]];
      q6Down = ArcTan[R36[[3, 1]], -R36[[3, 2]]];
      {{q4Up, q5Up, q6Up}, {q4Down, q5Down, q6Down}}], {joint, 
      jointAngles123}];
   (*Step 8:Combine all solutions*)
    allSolutions = 
    Flatten[Table[{Join[jointAngles123[[i]], q456Sol[[i, 1]]], 
       Join[jointAngles123[[i]], q456Sol[[i, 2]]]}, {i, 
       Length[jointAngles123]}], 1];
   (*Return all feasible joint angle configurations*) 
   allSolutions];

The result changes if I change the denavit hartenber parameters from
dhParameters = {
   {q1, d1, a1, -\[Pi]/2},
   {q2 - \[Pi]/2, 0, a2, 0},
   {q3, 0, a3, -\[Pi]/2},
   {q4, d4, 0, \[Pi]/2},
   {q5 + \[Pi], 0, 0, \[Pi]/2},
   {q6, d6, 0, 0}
   };
to 
dhParameters = {
   {q1, d1, a1, -\[Pi]/2},
   {q2 - \[Pi]/2, 0, a2, 0},
   {q3, 0, a3, -\[Pi]/2},
   {q4, d4, 0, \[Pi]/2},
   {q5 + \[Pi], 0, 0, -\[Pi]/2},
   {q6, d6, 0, 0}
   };