Matlab - ??? Index exceeds matrix dimensions.

Bueno, mi problema es el siguiente: tengo el código que os pongo más abajo.
Al meterle los parámetros según me va pidiendo con los datos del problema que tengo que resolver (Un flujo de carga según el método de Newton-Raphson) el programa se ejecuta sin problemas.
Ahora bien, yo obtengo al final del todo, cuatro vectores, que son Hd, Nd, Md y Ld.
Lo que quiero es crear cuatro matrices, cuya diagonal sea un vector de esos cuatro, es decir, una matriz con la diagonal que sea Hd, otra con la diagonal que sea Ns y así con los otros dos.
El problema viene cuando al ponerle:
diag(Hd) ó diag(Nd) ó diag(Md) ó diag(Ld), MatLab me muestra el error:

>> diag(Hd)
??? Index exceeds matrix dimensions.

No sé a qué es debido, y cómo solucionarlo y ya me estoy volviendo loco.
Podéis ayudarme??

Los datos que va pidiendo el programa (en este orden) son:
3

0
1
1
1
0
1
1
1
0

1
0.1j
0.2j
0.1j
1
0.25j
0.2j
0.25j
1

1
1
1

0
0
0

0
-1.5
1

0
-1.8
0

El código es el siguiente:

%%%%%%%%%%%%%%%CÓDIGO NEWTON-RAPHSON %%%%%%%%%%%%%%%%%

N = input ('Nº de nudos: ');

disp ('Matriz de relaciones entre nudos: 1 si hay relación, 0 si no hay relación ')
for i=1:N
for j=1:N
str=['Nudo ' num2str(i) ' y Nudo ' num2str(j) ': '];
R(i,j)=input(str);
end
end
R; %Matriz de relaciones "R". Son las relaciones que existen entre los nudos

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ESTO NO SIRVE PARA EL CÓDIGO, PERO SON COSAS MIAS
aa=0.02+0.04j;
bb=Inf;
cc=0.04+0.03j;
dd=0.01+0.01j;
ee=0.03+0.04j;
ff=0.1+0.1j;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp ('Matriz de impedancias; En la matriz de impedancias, las de las posiciones i=j se ponen igual a 1')
%Introducir la matriz de impedancias
%En la matriz de impedancias, las de las posiciones i=j se ponen igual a 1
%y donde no hay relación, se pone "Inf"
for i=1:N
for j=1:N
str=['Nudo ' num2str(i) ' y Nudo ' num2str(j) ': '];
Z(i,j)=input(str);
end
end
disp ('Matriz de impedancias "Z"')
Z; %Matriz de impedancias "Z"

m=1./Z;
y=m-eye(N);%Matriz de admitancias "y"

d=-y; %Matriz negativa de la de admitancias
v = [ ]; % vector vacío creado nuevo
for i = 1:N
v(i) = sum(y(i,:)); % Vector formado por la suma de las filas de la matriz de admitancias
end
diag=diag(v); %Matriz diagonal cuyos elementos son el vector "v"

disp ('Matriz de admitancias nodales "Y"')
Y=d+diag % Matriz de admitancias nodales "Y"
G=real(Y); %Matriz formada por la parte REAL de la matriz de admitancias nodales
B=imag(Y); %Matriz formada por la parte IMAGINARIA de la matriz de admitancias nodales

disp ('Tensiones de los nudos (con "perfil plano")')
for i=1:N
str=['Nudo ' num2str(i) ': '];
Vt(i)=input(str);
V=Vt'; %Vector de tensiones en los nudos en columna
end
V

disp ('Ángulos de los nudos (con "perfil plano")')
for i=1:N
str=['Tensión en el nudo ' num2str(i) ': '];
At(i)=input(str);
A=At'; %Vector de ángulos en los nudos en columna
end
A

disp ('Potencia activa (generada o cosumida) en los nudos (si es desconocida P = 0)')
for i=1:N
str=['P' num2str(i) ': '];
Pt(i)=input(str);
P=Pt'; %Vector de potencias activas en columna
end
P

disp ('Potencia reactiva consumida en los nudos (si es desconocida Q = 0)')
for i=1:N
str=['Q' num2str(i) ': '];
Qt(i)=input(str);
Q=Qt'; %Vector de potencias reactivas en columna
end
Q

% Tolerancia de precisión establecida en e=0.001
%e = input ('Tolerancia de precisión ');%Si se desea modificar en cada
%caso, sólo hay que eliminar el primer % de la lína 83 y borrar la línea 85
e=0.001;

m=G.*R; %Matriz para obtener los factores que intervienen en cada operación
n=B.*R; %Matriz para obtener los factores que intervienen en cada operación
[Aa,Ab]=meshgrid(A,A); %Matriz de los ángulos
Tt=Ab-Aa; %Matriz para calcular la resta de ángulos
T=(Tt)';

hp=m.*cos(T)+n.*sin(T); %Matriz que calcula la expresión Vj(Gij*Cos(Aij)+Bij*Sen(Aij))
hsump=(sum(hp))'; %vector suma de los valores de las columnas transpuesto
%%%%%disp ('Potencia activa calculada en los nudos conocidos "Pcal". Las incógnitas se muestran como "NaN". ')
Pcal=(V).*(V.*hsump);

hq=m.*sin(T)-n.*cos(T); %Matriz que calcula la expresión Vj(Gij*Sin(Aij)-Bij*Cos(Aij))
hsumq=(sum(hq))'; %vector suma de los valores de las columnas transpuesto
%%%%%%disp ('Potencia reactiva calculada en los nudos conocidos "Qcal". Las incógnitas se muestran como "NaN". ')
Qcal=(V).*(V.*hsumq);


%%%%%%disp ('Variación entre la potencia especificada (P,Q) y la potencia calculada (Pcal,Qcal)')
IncP=P-Pcal;
IncQ=Q-Qcal;
absIncP=abs(IncP); %Valor absoluto del Incremento de P
absIncQ=abs(IncQ); %Valor absoluto del Incremento de Q

%Cálculo del Jacobiano completo

%Diagonal de la matriz H
Hd=-Qcal-(((sum(B.*eye(size(B))))').*(V.^2));
%Diagonal de la matriz N
Nd=Pcal+(((sum(G.*eye(size(G))))').*(V.^2));
%Diagonal de la matriz M
Md=Pcal-(((sum(G.*eye(size(G))))').*(V.^2));
%Diagonal de la matriz L
Ld=Qcal-(((sum(B.*eye(size(B))))').*(V.^2));

Revisa que tus vectores tengan los valores y las dimensiones esperadas.

1 - Teclea Hd, a ver qué valores te da Matlab.
2 - También teclea size(Hd), para ver si realmente es vector, y no matriz o cell-array.
3 - Introduce los valores de Hd a mano, y luego teclea diag(Hd), a ver si funciona como lo esperas.

Hay varias maneras de depurar el código... pero revisa sección por sección, individualmente, no todo el programa completo, porque es demasiado extenso.

Es mejor separar modularmente el código, en porciones que se puedan probar individualmente. Cuando vayan funcionando bien, las incluyes en un código global que sabes que funciona bien parte por parte.

Saludos,

Voy a ir comprobando el código por partes a ver que saco.
Gracias por la ayuda y ya iré comentando que tal va la cosa.
Un saludo

Yo tenia un problema similar; el problema era que cuando se ingresan las matrices se tiene que dejar un espacio entre cada numero que ingreses por que de lo contrario el compilador te va a tomar como si fuera solo una cantidad.