This code calculates the load flow based on newton raphson methd for three bus power system. The Jacobian is written in a very easy form to understabd.

Newton Raphson Load Flow solving 3 Bus using MATLAB

%Find Ybus Admittance Matrix(any number of Bus)
%I am used 3 Bus and , 2 bus ans 3 bus are not connected
clc
clear
% |FromBus|ToBus|Impedance|LineCharging|
d = [ 1 2 .4j -20j ;
1 3 .25j -50j;
2 3 0j 0j ; ];
fb = d(:,1); tb = d(:,2);
y = 1./d(:,3);
% solving line Charginding don't used
if length(d)==3
b = y;
else
b = ((1./d(:,4))) + y;
end
nb = max(max(fb,tb));
Y = zeros(nb);
%put the value of diagonal position
Y(nb*(fb-1)+tb)=-y;
Y(nb*(tb-1)+fb)=-y;
%put the value diagonal position
for m=1:nb
Y(m,m)=sum(b(fb==m))+sum(b(tb==m));
end
%infinity problem solve (Inf or -Inf)
for m=1:nb*nb
if(real(Y(m))==Inf|real(Y(m))==-Inf)
Y(m)=0;
end
end

Newton Raphson Load Flow solving

% Newton Raphson Load Flow solving 3 Bus
clc
clear
%Finf out Ybus Admittance Matrix
% |FromBus|ToBus|Impedance|LineCharging|
d = [ 1 2 0.05j ];
fb = d(:,1); tb = d(:,2);
y = 1./d(:,3);
% solving line Charginding don't used
if length(d)==3
b = y;
else
b = ((1./d(:,4))) + y;
end
nb = max(max(fb,tb));
Y = zeros(nb);
%put the value of diagonal position
Y(nb*(fb-1)+tb)=-y;
Y(nb*(tb-1)+fb)=-y;
%put the value diagonal position
for m=1:nb
Y(m,m)=sum(b(fb==m))+sum(b(tb==m));
end
%infinity problem solve (Inf or -Inf)
for m=1:nb*nb
if(real(Y(m))==Inf|real(Y(m))==-Inf)
Y(m)=0;
end
end
Y_magnitude =abs(Y);

Y_theta=angle(Y)
%gavien valu
v1=1.0;theta1=0;pq2=1+0.5j;
p2=-real(pq2);q2=-imag(pq2);
%inatial guess
v2=1;theta2=0;
X=[v2 theta2]';
%
for i=1:200
fp2=Y_magnitude(2,1)*v1*v2*cos(Y_theta(2,1)+theta1-theta2)-p2;
fq2=-Y_magnitude(2,1)*v1*v2*sin(Y_theta(2,1)+theta1-
theta2)+Y_magnitude(2,2)*v2^2-q2;
%jacobian matrix find
j(1,1)=Y_magnitude(2,1)*v1*v2*sin(Y_theta(2,1)+theta1-theta2);
j(1,2)=Y_magnitude(2,1)*v1*cos(Y_theta(2,1)+theta1-theta2);
j(2,1)=Y_magnitude(2,1)*v1*v2*cos(Y_theta(2,1)+theta1-theta2);
j(2,2)=-Y_magnitude(2,1)*v1*sin(Y_theta(2,1)+theta1-
theta2)+2*Y_magnitude(2,2)*v2;
%Newton Raphson applied
X=X-inv(j)*[fq2 fp2]';
theta(i)=X(2);
voltage(i)=X(1);
theta2=X(2);v2=X(1);
end
v2
theta2