matriz admitancias

Exercício que apresenta os cálculos de matriz de impedâncias e potência instantânea

Nota: Todos os valores das impedâncias estão em Ω.
É preciso reparar que a impedância Z2 é definida como o sumatório das duas impedâncias (a resistiva e a capacitiva), e não como o sumatório dos inversos das impedâncias.
Também que a impedância Z5 e a impedância Z6 são independentes. Isto quer dizer que no nó onde se encontram as impedâncias Z2, Z5 e Z6, o valor do sumatório das admitâncias é:
Z2 = 2+3j Ω
Y_22 = 1/Z2 + 1/Z5 + 1/Z6 + 1/Z3
Consequentemente, no nó onde se encontram as impedâncias Z5, Z6, Z4 e Z3, o valor do sumatório das admitâncias é:
Y_33 = 1/Z5 + 1/Z6 + 1/Z4
A expressão da potência instantânea: p = VIr (1+ cos ( 2 (ωt + δ) ) ) + V*Ix * sen ( 2* ( ωt + δ ) )
O programa em Matlab ou Octave que resolve o problema está a seguir:
% Aula 2 SEE Laboratorio
%Definicoes de variaveis
Z1 = 3 - 10*i;
Z2 = i*3+2;

Z3 = -i+2;
Z4 = 2-4*i;
Z5 = -2*i;
Z6 = i;

%Definicao componentes matriz de admitancias
Y_11 = 1/Z1 + 1/Z2;
Y_22 = 1/Z2 + 1/Z5 + 1/Z6 + 1/Z3;
Y_33 = 1/Z5 + 1/Z6 + 1/Z4;
Y_12 = -1/Z2;
Y_21 = Y_12;
Y_23 = -(1/Z5 + 1/Z6);
Y_32 = Y_23;
Y_13 = 0;
Y_31 = Y_13;
%Matriz de admitancias
Y = [Y_11 Y_12 Y_13;Y_21 Y_22 Y_23;Y_31 Y_32 Y_33];

%Definicao componentes vector de correntes
I1 = 0.5*exp(-i*(10*pi/180));
I2 = 0.8*exp(-i*(35*pi/180));
I3 = 0.7*exp(i*(10*pi/180));
%Vector de correntes
I = [I1;I2;I3];
%Solucao de tensoes para os nos
V = Y^-1*I;
%Corrente na impedancia Z4
I4 = V(3)/Z4;
%Potencia aparente na impedancia Z4
S4 = V(3)*conj(I4);
%Desenho potencia consomida pela impedancia Z4
%Assumindo uma frequencia de 50Hz
%Fase da tensao

delta = angle(V(3));
%Fase da corrente beta = angle(I4);
I_r = norm(I4)*cos(delta-beta);
I_x = norm(I4)*sin(delta-beta); freq = 50;
T_max = 1/freq; omega_ = 2*pi*freq; time=0:0.0001:T_max/2; % p = V*I_r (1+ cos ( 2* (omega*time + delta) ) ) + V*I_x * sen ( 2* (omega*time + delta