潮流计算之前推回代法的matlab程序

function f=myfile（nod8） 
clc
clear all;
format long;
Z=[0  1  0.54+j*0.65  0.6+j*0.45
   1  2  0.62+j*0.5  0.4+j*0.3
   2  3  0.6+j*0.35  0.4+j*0.28
   1  4  0.72+j*0.75  0.6+j*0.4
   4  5  1+j*0.55  0.4+j*0.3
   4  6  0.65+j*0.35  0.5+j*0.35
   1  7  0.9+j*0.5    0.5+j*0.4];%33nodes datas 第一列为首节点编号 第二列为末节点编号 第三列为支路参数 第四列为节点负荷

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%开始计算
V=ones（71）;%赋节点电压值
for i=1:7
    I=（Z（:4）.‘）*（diag（1./V））;
    %V1=V;
    I=I‘;
    I1=I;
   for k=7:-1:1   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%后推支路电流及功率分布 
      A=（find（Z（:1）==k））;%判断矩阵Z第1列是否存在值为k的元素，并把行值赋给矩阵A，亦即找出网络的末节点
        if size（A）==0  %为末节点时
           I（k1）;      
        else  %不为末节点
           I（k1）=sum（I（A1））+I（k1）; %利用KCL求支路电流
        end
   end
   V（11）=10.5-Z（13）*I（11）;%%%%%%%%%%%