extending matrix with numbers (not zero)

In the folowing code; I have error

Index in position 1 exceeds array bounds (must not exceed 4).

Error in Untitled1111111 (line 137)

load(1,(head(1,i)-1)*(dof+j))=fload(j+1,i);

the matrices dimention is:

load(1x300), head (4x30), fload(4x30)

and the cod is:

close all;
clear all
dbstop if error
% Initial data:
TubeLength =2.38e-9;
%%TubeLength = input ('TubeLength'); %m 
TubeDiameter=1.62e-9;
%%%TubeDiameter =  input ('TubeDiameter'); %m
r=TubeDiameter/2;
A=pi*r^2;
D = 3; %degrees of freedom at each node
presc_disp=1e-9; %m
L=0.141e-9; %m Length of a C-C bond
numIterations=1; %number of iterations
numNaN=0;
% %%%%%  the ref .
%%Force constants from Tserpes et al. (2005):
% Kr is EA=6.52e-7 %N·nm^-1
% ktheta is EI=8.76e-10 %N·nm·rad^-2
% ktao is GJ=2.78e-10 %N·nm·rad^-2
EA=L*652 %N
EI=L*8.76e-19 %N·m^2·rad^-2
GJ=L*2.78e-19; %N·m^2·rad^-2
% Transformation to nondimensional space:
% [F]=EI/L^2; [L]=L. Then, multiply the resulting forces with
% EI/L^2 and displacements with L to get forces and
% displacements in nm.
Lreal=L;
EIreal=EI;
presc_disp_real=presc_disp;
EA=EA/EI*L^2;
GJ=GJ/EI;
presc_disp=presc_disp/L;
L=1;
EI=1;
% Stiffness matrix of a C-C bond (adapted to 3 degree of freedom):
KBond=[EA/L 0 0 -EA/L 0 0 
0 12*EI/L^3 0 0 -12*EI/L^3 0 ;
0 0 12*EI/L^3 0 0 -12*EI/L^3;
-EA/L 0 0 EA/L 0 0 ;
0 -12*EI/L^3 0 0 12*EI/L^3 0 ;
0 0 -12*EI/L^3 0 0 12*EI/L^3]
% Calculations:


[elementss1]=xlsread('elementss1.xlsx');
[nodess1]=xlsread('nodess1.xlsx');
%%nnodess = input ('number  of nodes');
%%nElementss = input ('number of elements');
nnodess1= 200;               %%% number of nodes
nElementss2 = 200;              %%%% number of elements
nodesDef=zeros(1,nnodess1);

k=zeros(nnodess1*D,nnodess1*D);
% Defects:
defectsPercentage=1; %per cent
numDef=round((nnodess1*defectsPercentage)/100);
for i=1:numDef
    def=round((nodess1-1)*rand)+1; %The addition of 1 and -1 is to avoid zero as the resulting random number. elements(D,:)=[];
nnodess1=nnodess1-1;
end
zeroMatrix = [0 0 0; 0 0 0; 0 0 0];
vi = [1 0 0];
vj = [0 1 0];
vk = [0 0 1];
for i=1:nnodess1-1
n = nodess1(2,i);
m = nodess1(3,i)
localZ = [1-2, 3-2, 3-4];
localZ = localZ/sqrt(localZ*localZ');
% Find local y
if not(localZ(3)==0)
y3 = (-localZ(1)-2*localZ(2))/localZ(3);
localY = [1 2 y3];
elseif not(localZ(2)==0)
y2 = (-localZ(1)-2*localZ(3))/localZ(2);
localY = [1 y2 2];
elseif not(localZ(1)==0)
y1 = (-localZ(2)-2*localZ(3))/localZ(1);
localY = [y1 1 2];
end
localY = localY/sqrt(localY*localY');
% Find local z
localX = cross(localY, localZ);
localX = localX/sqrt(localX*localX');
% Values for rotation matrix
cosZI = localZ*vi';      
cosZJ = localZ*vj';
cosZK = localZ*vk';
cosYI = localY*vi';
cosYJ = localY*vj';
cosYK = localY*vk';
cosXI = localX*vi';
cosXJ = localX*vj';
cosXK = localX*vk';
rotationSmall = [
cosZI cosZJ cosZK;
cosYI cosYJ cosYK;
cosXI cosXJ cosXK];
RR = [rotationSmall zeroMatrix;zeroMatrix rotationSmall];
KElement = abs(RR'*KBond*RR)

n = nodess1(2,i);
m = nodess1(3,i);
k((n-1)*D+1:n*D,(n-1)*D+1:n*D)=k((n-1)*D+1:n*D, (n-1)*D+1:n*D)+KElement(1:D, 1:D);
k((m-1)*D+1:m*D,(m-1)*D+1:m*D)=k((m-1)*D+1:m*D, (m-1)*D+1:m*D)+KElement(D+1:2*D, D+1:2*D);
k((n-1)*D+1:n*D,(m-1)*D+1:m*D)=k((n-1)*D+1:n*D, (m-1)*D+1:m*D)+KElement(1:D, D+1:D*2);
k((m-1)*D+1:m*D,(n-1)*D+1:n*D)=k((m-1)*D+1:m*D, (n-1)*D+1:n*D)+KElement(D+1:D*2, 1:D);

end
head = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494 500 501 503 504;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
];
fload = [1 4 5 16 17 34 35 51 52 72 73 86 87 476 477 478 479 482 483 488 489 490 491 492 493 494 500 501 503 504;
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 ;
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 ;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
]*presc_disp;
dof=3;nNodes=100;
save k_real.mat k;
load=zeros(1,dof*nNodes);
for i=1:size(head(1,:),2)
for j=1:dof;
if head(j+1,i)~=0 % displacement prescribed
 
    
for m=1:dof*nNodes;
load(m)=load(m)-k(m,(head(1,i)-1)*dof+j)*fload(j+1,i);
end
i=i+1;
j=j+1;
k((head(1,i))*(dof+j),:)=zeros(1,dof*2*nNodes);
k(:,(head(1,i)-1)*dof+j)=zeros(dof*2*nNodes,1);
k((head(1,i)-1)*dof+j,(head(1,i)-1)*dof+j)=1.0;
load(1,(head(1,i)-1)*(dof+j))=fload(j+1,i);
if fload(j+1,i) ~= 0
kdof=j;



end
end
end
end


clear RR KEelement;
clear elements;
disps=loads/k;
clear k;
loads k_real.mat
forces=k*disps'
clear k
tip=[476 477 478 479 482 483 488 489 490 491 492 493 494];
% 'tip'contains the node numbers at one tip
avdisp=0;
for i=1:size(tip,2)
tip_disps(i)=disps((tip(i)-1)*dof+kdof);
avdisp=avdisp+disps((tip(i)-1)*dof+kdof);
end
fixed=[1 4 5 16 17 34 35 51 52 72 73 86 87];
% 'fixed'contains the node numbers at the other tip
totforce=0;
for i=1:size(tip,2)
totforce=totforce+forces((tip(i)-1)*dof+kdof);
end
avdisp=avdisp/size(tip,2);
avdisp=avdisp*Lreal; % in meters
tip_disps=tip_disps*Lreal; % in meters
totforce = totforce * EIreal/Lreal^2; % in Newtons
YY=totforce/(TubeDiameter*avdisp/TubeLength*pi) % N/m=Pa*m
Young(it)=YY;
it=it+1;