//Determine the voltages starting with a flat voltage profile.
clear 
clc;

Y=[3-%i*12 -2+%i*8 -1+%i*4 0;-2+%i*8 3.666-%i*14.664 -.666+%i*2.6664 -1+%i*4;-1+%i*4 -.666+%i*2.6664 3.666-%i*14.664 -2+%i*8;0 -1+%i*4 -2+%i*8 3-%i*12];
P2=.5;
P3=-.4;
P4=-.3;
Q4=-.1;
Q3=-.3;
V3=1;
V4=1;
V1=1.06;
V2=1.04;
V30=1;
V40=1;
Q2=-imag([V2*[Y(2,1)*V1+Y(2,2)*V2+Y(2,3)*V3+Y(2,4)*V4]]);
V21=(((P2-%i*Q2)/V2)-Y(2,1)*V1-Y(2,3)*V30-Y(2,4)*V40)/(Y(2,2));
d=atand(0.0291473/1.0472868);
V21=1.04*(cosd(d)+%i*sind(d));
disp(V21,"V21=");
V31=(((P3-%i*Q3)/V3)-Y(3,1)*V1-Y(3,2)*V21-Y(3,4)*V40)/(Y(3,3));
disp(V31,"V31=");
V41=(((P4-%i*Q4)/V4)-Y(4,2)*V21-Y(4,3)*V31)/(Y(4,4));
disp(V41,"V41=");