I have checked everything again and I'm pretty confident that I didn't make a mistake. The cause of the difference is actual theta offset i.e. what is considered a variable. For example, when setting up coordinate systems, I have considered that zero positions of all angles corresponds to the same configuration that you considered as: [0, -90, 90, 0, 0, 0].
So when you calculate transformation matrix for [0,-60,90,0,0,0], I need to calculate for [0, +30, 0, 0, 0, 0], because comparing to mine assumed zero configuration, only th2 is changed for +30° because Z axis of joint 2 is oriented "into the screen from us".
So the following code:
% D_H parameters: A1d = [0, 30, 0, 0.0, 0.0, 0.0]; th = deg2rad(A1d); a1 = 500; alpha1 = -pi/2; d1 = 1045; th1 = th(1); a2 = 1300; alpha2 = 0; d2 = 0; th2 = th(2)-pi/2; a3 = 55; alpha3 = -pi/2; d3 = 0; th3 = th(3)+pi; a4 = 0; alpha4 = pi/2; d4 = -1025; th4 = th(4); a5 = 0; alpha5 = -pi/2; d5 = 0; th5 = th(5); a6 = 0; alpha6 = pi; d6 = -290; th6 = th(6); T1 = DH_calc(a1,alpha1,d1,th1); T2 = DH_calc(a2,alpha2,d2,th2); T3 = DH_calc(a3,alpha3,d3,th3); T4 = DH_calc(a4,alpha4,d4,th4); T5 = DH_calc(a5,alpha5,d5,th5); T6 = DH_calc(a6,alpha6,d6,th6); T06_dh = T1*T2*T3*T4*T5*T6; vpa(T06_dh,6)
will produce the following output: