I have created a newton raphson algorithm based load flow programming for 30 bus system.After run it shows unexpected matlab expression but it does not show the line no where the fault has occurred. How can I solve it?

zdata = [1 30 0.0192 0.0575 0.02640 1 1 3 0.0452 0.1852 0.02040 1 30 4 0.0570 0.1737 0.01840 1 3 4 0.0132 0.0379 0.00420 1 30 29 0.0472 0.1983 0.02090 1 30 6 0.0581 0.1763 0.01870 1 4 6 0.0119 0.0414 0.00450 1 29 7 0.0460 0.1160 0.01020 1 6 7 0.0267 0.0820 0.00850 1 6 28 0.0120 0.0420 0.00450 1 6 9 0.0 0.2080 0.0 0.978 6 10 0 0.5560 0 0.969 9 27 0 0.2080 0 1 9 10 0 0.1100 0 1 4 12 0 0.2560 0 0.932 12 26 0 0.1400 0 1 12 14 0.1231 0.2559 0 1 12 15 0.0662 0.1304 0 1 12 16 0.0945 0.1987 0 1 14 15 0.2210 0.1997 0 1 16 17 0.0824 0.1923 0 1 15 18 0.1073 0.2185 0 1 18 19 0.0639 0.1292 0 1 19 20 0.0340 0.0680 0 1 10 20 0.0936 0.2090 0 1 10 17 0.0324 0.0845 0 1 10 21 0.0348 0.0749 0 1 10 22 0.0727 0.1499 0 1 21 22 0.0116 0.0236 0 1 15 23 0.1000 0.2020 0 1 22 24 0.1150 0.1790 0 1 23 24 0.1320 0.2700 0 1 24 25 0.1885 0.3292 0 1 25 13 0.2544 0.3800 0 1 25 11 0.1093 0.2087 0 1 8 11 0 0.3960 0 0.968 11 5 0.2198 0.4153 0 1 11 2 0.3202 0.6027 0 1 5 2 0.2399 0.4533 0 1 28 8 0.0636 0.2000 0.0214 1 6 8 0.0169 0.0599 0.065 1]; Y=fybus(zdata); V=[1.06;1;1;1.06;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1.071;1.082;1.01;1.01;1.043]; d=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]; Ps=[-10.6;-2.4;-7.6;-2.4;-0.0;-22.8;-0.0;-0.0;-5.8;-0.0;-11.2;-3.5;-6.2;-8.2;-3.5;-9.0;-3.2;-9.5;-2.2;-17.5;-0.0;-3.2;-8.7;-0.0;-0.0;-0.0;-30.0;-94.2;18.3]; Qs=[-1.9;-1.2;-1.6;-0.9;-0.0;-10.9;-0.0;-0.0;-2.0;-0.0;-7.5;-2.3;-1.6;-2.5;-1.8;-5.8;-0.9;-3.4;-0.7;-11.2;-0.0;-1.6;-6.7;-0.0;-0.0;-0.0;-30.0;-19.0;-12.7]; y = abs(Y); t = angle(Y); iter=0; pwracur = 0.00025; DC=10; N=30; M=3; while max(abs(DC)) > pwracur; iter = iter + 1; P=zeros(N-1,1); k=0; for i=2:N; k=k+1; P(k)=0;

    for j=1:N;

       Pv=V(i)*V(j)*y(i,j)*cos(t(i,j)-d(i)+d(j));

       P(k)=P(k)+Pv; 

    end
 end
P;
Q=zeros(N-M-1,1);
k=0;
for i=[2:N-M];
      k=k+1;
   Q(k)=0;
     for j=[1:N];
           Qv=-V(i)*V(j)*y(i,j)*sin(t(i,j)-d(i)+d(j));
           Q(k) = Q(k)+Qv;
      end
end
Q;
J1=zeros(N-1,N-1);
for i=2;
   for k=1:N-1;
        J1(k,k)=0;
      for j=1:N;
        if i~=j;
        J1(k,k)=J1(k,k)+V(i)*V(j)*y(i,j)*sin(t(i,j)-d(i)+d(j));
        else
        end
      end
    a=a+1;
  end 
end
for i=2;
  for k=1:N-1;
    for j=2:N;
      for n=1:N-1;
        if i~=j && j==n+1;
          if k~=n;
        J1(k,n)=-V(i)*V(j)*y(i,j)*sin(t(i,j)-d(i)+d(j));
       else,end
       else,end
      end
    end
    i=i+1;
  end
end
J2=zeros(N-1,N-1-M);
for i=2;
   for k=1:N-1-M;
         J2(k,k)=2*V(i)*y(i,i)*cos(t(i,i));
      for j=1:N;
       if i~=j;
        J2(k,k)=J2(k,k)+V(j)*y(i,j)*cos(t(i,j)-d(i)+d(j));
      else,end
     end
    i=i+1;
   end
end
for i=2;
    for k=1:N-1;
      for j=2:N;
        for n=1:N-1-M;
            if i~=j;
            if k~=n;             
            J2(k,n)=V(i)*y(i,j)*cos(t(i,j)-d(i)+d(j));
        else,end
        else,end
       end
      end
     i=i+1;
   end
end
J3=zeros(N-1-M,N-1);
for i=2;
   for k=1:N-1-M;
    J3(k,k)=0;
     for j=1:N;
        if i~=j;
         J3(k,k)=J3(k,k)+V(i)*V(j)*y(i,j)*cos(t(i,j)-d(i)+d(j));
        else,end
      end
    i=i+1;
   end
end
for i=2;
    for k=1:N-1-M;
      for j=2:N;
        for n=1:N-1;
         if i~=j & j==n+1;
           if k~=n;
         J3(k,n)=-V(i)*V(j)*y(i,j)*cos(t(i,j)-d(i)+d(j));
        else,end
        else,end
      end
      end
    i=i+1;
  end
end
J4=zeros((N-1-M),(N-1-M)); 
for i=2;
    for k=1:N-1-M;
         J4(k,k)=-2*V(i)*y(i,i)*sin(t(i,i));
       for j=1:N;
         if i~=j;
         J4(k,k)=J4(k,k)-V(j)*y(i,j)*sin(t(i,j)-d(i)+d(j));
        else,end
      end
    i=i+1;
  end
end
for i=2;
 for k=1:N-1-M;
   for j=2:N;
     for n=1:N-1-M;
       if i~=j;
          if k~=n;
         J4(k,n)=-V(i)*y(i,j)*sin(t(i,j)-d(i)+d(j));
       else,end
       else,end
      end
   end
    i=i+1;
  end
end
J =[J1 J2;J3 J4];
DP = Ps - P;
DQ = Qs - Q; 
DC = [DP ; DQ];
DX = J\DC;
for i=2:N;
 d(i)= d(i)+DX(i-1);
end
L=length(d);
for i=2:N-M;
  V(i)=V(i)+DX(i+L-2);
end
delta =180/pi*d;
Result=[V, d, delta];
end
for i=1:N;Q(i)=0;
    for j=1:N
Q(i)= Q(i)-V(i)*V(j)*y(i,j)*sin(t(i,j)-d(i)+d(j));
end
end
for j=1:N
   P(j)=0;
Q(j)=0;
for i=1:N; 
P(j)= P(j)+V(j)*V(i)*y(j,i)*cos(t(j,i)-d(j)+d(i));
Q(j)=Q(j)-V(j)*V(i)*y(j,i)*sin(t(j,i)-d(j)+d(i));
end
end
P; 
Q;
disp('   Voltage    angle     P          Q');
fprintf('%g'),disp([            V                d                P                    Q])