Example code Load Flow using fsolve Matlab

clear
clc

V1=1.01;

ZL12=0.02+i*0.01;
ZL23=0.03+i*0.01;
Zload2=9.8+i*1.96;
Sload3=0.2+i*0.05;

YL12=1/ZL12;
YL23=1/ZL23;
Yload2=1/Zload2;
Pload3=real(Sload3);
Qload3=imag(Sload3);

Y=[YL12,-YL12,0;
-YL12,YL12+YL23+Yload2,-YL23;
0,-YL23,YL23];

G12=real(Y(1,2));
B12=imag(Y(1,2));

G22=real(Y(2,2));
B22=imag(Y(2,2));

G23=real(Y(2,3));
B23=imag(Y(2,3));

G33=real(Y(3,3));
B33=imag(Y(3,3));

f=@(V2,V3,theta2,theta3) [V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3)); V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+V3*V3*G33+Pload3; V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3)); V2*V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-V3*V3*B33+Qload3];

fp=@(x) f(x(1),x(2),x(3),x(4));

[x, fval, info] = fsolve (fp, [1.01;1.01;0;0]);

disp(x);

V2=x(1);
V3=x(2);
theta2=x(3);
theta3=x(4);