SOlving power flow using gauss seidel and matlab, HELP

[O.J.]

Homework Statement

We were asked to do a course project using to solve a powerflow system usign gauss seidel and Newton raphson methods. after writing down the equations at each bus and double checking them this is the final code2. The attempt at a solution
c=1
V1=1
P2=-8
Q2=-2.8
V3=1.05
P3=4.4
P4=0
Q4=0
P5=0
Q5=0
V2=1
V4=1
V5=1
th2=atan (imag(V2)./real(V2))
th3=atan (imag(V3)./real(V3))
th4=atan (imag(V4)./real(V4))
th5=atan (imag(V5)./real(V5))
error = 0.001
errorm = 1

while (errorm > error)

V2new = (1./(2.678-i*28.459))*((8+i*2.8)-((-0.893+i*0.919)*V4+(-1.785+i*19.839)*V5))
error2 = abs(V2new - V2)th3new = acos(real(((11.92*sin(th4)-147.96*cos(th4)) + V2*((-0.89*sin(th4-th2)-9.92*cos(th4-th2))+1.05*(cos (th3)+i*sin (th3))*(-7.46*sin(th4-th3)-99.44*cos(th4-th3))+V4*147.96+V5*(-7.48*sin(th4-th5)-99.44*cos(th4-th5))))/4.19))
error3 = abs(th3new - th3)V4new = (1./(11.92-i*147.96)) * (-1*((-0.89+i*9.92)*V2+(-7.46+i*99.44)*1.05*(cos(th4)+i*sin(th3))+(-3.57+i*39.68)*V5))
error4 = abs(V4new - V4)V5new = (1./(9.09-i*107.86)) * (-1*((-3.74+i*49.72)+V2*(-1.786+i*19.839)++V4*(-3.72+i*39.68)))
error5 = abs(V5new - V5)

V2 = V2new
th3 = th3new
V4 = V4new
V5 = V5new
th2=atan (imag(V2)./real(V2))
th4=atan (imag(V4)./real(V4))
th5=atan (imag(V5)./real(V5))

error_1 = max(error2, error3)
error_2 = max(error4, error5)
errorm = max(error_1, error_2)
c=c+1
end
this is my code. please try it and tell me anything if you notice anything.

 

[O.J.]

ok, I am tracing my output to try to know where my mistakes are. i discovered that my expresion for the error was a mistake, because it only took into account the real part so my new expression is c=1
V1=1
P2=-8
Q2=-2.8
V3=1.05
P3=4.4
P4=0
Q4=0
P5=0
Q5=0
V2=1
V4=1
V5=1
th2=atan (imag(V2)./real(V2))
th3=atan (imag(V3)./real(V3))
th4=atan (imag(V4)./real(V4))
th5=atan (imag(V5)./real(V5))
error = 0.0001
errorm = 1

while (errorm > error)

V2new = (1./(2.678-i*28.459))*((8+i*2.8)-((-0.893+i*0.919)*V4+(-1.785+i*19.839)*V5))
error2 = abs(((real(V2new).^2+imag(V2new).^2)^0.5)-((real(V2).^2+imag(V2).^2).^0.5))th3new = asin(imag(((11.92*sin(th4)-147.96*cos(th4)) + V2*((-0.89*sin(th4-th2)-9.92*cos(th4-th2))+1.05*(cos (th3)+i*sin (th3))*(-7.46*sin(th4-th3)-99.44*cos(th4-th3))+V4*147.96+V5*(-7.48*sin(th4-th5)-99.44*cos(th4-th5))))./4.19))
error3 = abs(th3new-th3)

V4new = (1./(11.92-i*147.96)) * (-1*((-0.89+i*9.92)*V2+(-7.46+i*99.44)*1.05*(cos(th4)+i*sin(th3))+(-3.57+i*39.68)*V5))
error4 = abs(((real(V4new).^2+imag(V4new).^2).^0.5)-((real(V4).^2+imag(V4).^2).^0.5))V5new = (1./(9.09-i*107.86)) * (-1*((-3.74+i*49.72)+V2*(-1.786+i*19.839)+V4*(-3.72+i*39.68)))
error5 = abs(((real(V5new).^2+imag(V5new).^2).^0.5)-((real(V5).^2+imag(V5).^2).^0.5))

V2 = V2new
th3 = th3new
V4 = V4new
V5 = V5new
th2=atan (imag(V2)./real(V2))
th4=atan (imag(V4)./real(V4))
th5=atan (imag(V5)./real(V5))

error_1 = max(error2, error3)
error_2 = max(error4, error5)
errorm = max(error_1, error_2)
c=c+1
end

Im still getting NaN results as errors are growing huge, I traced the output and found where I could be going wrong but I didnt resolve it so far. th3 is coming out as a complex number when it shuold be an angle only. I don't know wat I am doing wrong.. please help