Example code Load Flow using Newton Raphson 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];

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

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

x=zeros(3,4);
x(1,:)=[1.01,1.01,0,0];
disp(norm(fp(x(1,:))));
for n=2:length(x)
    x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';
    disp(norm(fp(x(n,:))));
end
disp('Solution:');
disp(x(n,:));