Someone help me debug this Matlab Functions initialisation code. It gives errors such as 'Inner matrix dimensions must agree' and so on. Any assistance will be gladly appreciated.

if true
% file name: initialization.m
% started by callback of Simulink model
% initialises data for buses grid,
% creates snet, admittance matrix,
% 50/60 Hz initialization (deltime)steptime
% calculates a first voltage bus estimation
j=sqrt(-1);
%%Set generation
GeneratorBus = [
% nr ... Vmag P Q ...
% 1, 0, 1.000, 0.5, 0.000, 0.000;
2, 0, 1.045, 0.483, 0.000, 0.000;
3, 0, 1.010, 0.342, 0.000, 0.000;
% 6, 0, 1.070, 0.112, 0.000, 0.000;
6, 0, 1.070, 0.2, 0.000, 0.000;
8, 0, 1.090, 0.200, 0.000, 0.000; % was zero Pm ouput before!
% 14, 0, 1.036, -0.08, 0, 0;
];
%%Set load
LoadBus = [
% nr ... P Q
4, 0, -0.4, -0.5;
5, 0, -0.07, -0.01;
7, 0, 0.000, 0.000;
9, 0, -0.35, -0.05; % was Q=-.15
10, 0, -0.20, -0.06;
11, 0, -0.2, -0.08;
12, 0, -0.16, -0.05; % was Q=-.08
13, 0, -0.28, -0.05; % was Q=-.08
14, 0, -0.08, 0.05; % was Q=-.05
];
%%Set swing, slack bus
SwingBus = [
% nr ... Vmag P Q ...
1, 0, 1.060, 0.000;
];
%%setlines
LineData = [
% from to impedance
2, 1, 4.999-15.263j;
5, 1, 1.026-4.235j;
3, 2, 1.135-4.782j;
4, 2, 1.686-5.116j;
5, 2, 1.701-5.194j;
4, 3, 1.986-5.069j;
5, 4, 6.841-21.579j;
7, 4, 0.000-4.782j;
9, 4, 0.000-1.798j;
6, 5, 0.000-3.968j;
11, 6, 1.955-4.094j;
12, 6, 1.526-3.176j;
13, 6, 3.099-6.103j;
8, 7, 0.000-5.677j;
9, 7, 0.000-9.090j;
10, 9, 3.902-10.365j;
14, 9, 1.424-3.029j;
11, 10, 1.881-4.403j;
13, 12, 2.489-2.252j;
14, 13, 1.137-2.315j;
];
LineMatrix = [
% from to
% to from
4.999-15.237*j, -4.999+15.263*j;
-4.999+15.263*j, 4.999-15.237*j;
1.026-4.210*j, -1.026+4.235*j;
-1.026+4.235*j, 1.026-4.210*j;
1.135-4.760*j, -1.135+4.782*j;
-1.135+4.782*j, 1.135-4.760*j;
1.686-5.097*j, -1.686+5.116*j;
-1.686+5.116*j, 1.686-5.097*j;
1.701-5.177*j, -1.701+5.194*j;
-1.701+5.194*j, 1.701-5.177*j;
1.986-5.052*j, -1.986+5.069*j;
-1.986+5.069*j, 1.986-5.052*j;
6.841-21.572*j, -6.841+21.579*j;
-6.841+21.579*j, 6.841-21.572*j;
0.000-4.782*j, 0.000+4.782*j;
0.000+4.782*j, 0.000-4.782*j;
0.000-1.798*j, 0.000+1.798*j;
0.000+1.798*j, 0.000-1.798*j;
0.000-3.968*j, 0.000+3.968*j;
0.000+3.968*j, 0.000-3.968*j;
1.955-4.094*j, -1.955+4.094*j;
-1.955+4.094*j, 1.955-4.094*j;
1.526-3.176*j, -1.526+3.176*j;
-1.526+3.176*j, 1.526-3.176*j;
3.099-6.103*j, -3.099+6.103*j;
-3.099+6.103*j, 3.099-6.103*j;
0.000-5.677*j, 0.000+5.677*j;
0.000+5.677*j, 0.000-5.677*j;
0.000-9.090*j, 0.000+9.090*j;
0.000+9.090*j, 0.000-9.090*j;
3.902-10.365*j, -3.902+10.365*j;
-3.902+10.365*j, 3.902-10.365*j;
1.424-3.029*j, -1.424+3.029*j;
-1.424+3.029*j, 1.424-3.029*j;
1.881-4.403*j, -1.881+4.403*j;
-1.881+4.403*j, 1.881-4.403*j;
2.489-2.252*j, -2.489+2.252*j;
-2.489+2.252*j, 2.489-2.252*j;
1.137-2.315*j, -1.137+2.315*j;
-1.137+2.315*j, 1.137-2.315*j;];
%%Admittanc matrix and snet vector
slacklist=SwingBus(:,1); %bus number which is the swing
loadlist=LoadBus(:,1); %bus numbers with loads
genlist=GeneratorBus(:,1); %similar...
nslack=length(slacklist);
ngen=length(genlist);
nload=length(loadlist);
%nbus=nslack+nload+ngen;
nbus = nload+ngen;
vbus=ones(nbus,1);
vbus(genlist)=GeneratorBus(:,3); %setting the voltages of gens
vbus
size(vbus)
size(yb)
%vbus(slacklist)=SwingBus(:,3);
snet=(1+j)*ones(nbus,1); %initialization with 1+j
snet(loadlist)=-LoadBus(:,3)-j*LoadBus(:,4);
%inverse of minus sign (load in now positive)
snet(genlist)=-GeneratorBus(:,4)-j*GeneratorBus(:,5);
%inverse (gen is now negative)
fromto=LineData(:,1:2);
index=(1:length(fromto));
nline=length(index); %equivalent to length(fromto)
index;
nline;
oddindex=(2*index)-(ones(1,nline));
evenindex=2*index;
new=[fromto(:,1);fromto(:,2)]; %column vector
new(oddindex,1)=fromto(:,1);
new(evenindex,1)=fromto(:,2);
nbus=max(new);
aa=zeros(nbus,nline);
for kk = 1:length(new)
aa(new(kk),kk)=1;
end
kkl=1:length(new);
aa=sparse(new,kkl,ones(length(new),1));
lkl=1:nline;
lkl2=1:(2*nline);
lkodd=2*lkl-ones(1,nline);
lkeven=2*lkl;
lktotal1=[lkodd lkeven lkodd lkeven];
lktotal2=[lkodd lkeven lkeven lkodd];
bigline=[LineMatrix(lkodd,1); LineMatrix(lkeven,2); ...
LineMatrix(lkodd,2);LineMatrix(lkeven,1);];
bb=sparse(lktotal1,lktotal2,bigline,40,40);
% "primitive" (element-by-element) admittance matrix for
% network, with each element treated as a two port
% (hence, for 20 transmission line elements, we
% get a block diagonal, 40x40 matrix. Each
% transmission line yields a 2x2 block along its diagonal.
% This somewhat non-standard formulation accomodates
% future additions of elements such as phase shifting
% transformers, which can not be represented as simple
% two terminal circuit elements.
yb=aa*bb.*aa; % yb is the bus admiitance matrix for the study network.
% Note that the incidence matrix, aa, and
% primitive admittance matrix, bb, are individualy
% stored. This allows the transient stability
% simulation routine ("gotrans.m") an easy
% mechanism for switching transmission lines
% in and out of the network.
yb=full(yb);
%%load data and modify data from original data:
snet(1)=-0.4951 +j*0; % slack node
snet(2)=-0.2; % 20090630
snet(3)=-0.4;
snet(6)=-0.5;
snet(8)=-0.2; % 20090630
%
%%60 hz 50 hz initialization
%deltime=.0166666; %60Hz
deltime=.02; %50Hz
%%constants for Swing equation (D (damping) and M (generator inertias))
% original: D=.004*ones(length(genlist),1); % generator rotational damping
syms D1 D2 D3 D6 D8 D14;
D=[D1; D2; D3; D6; D8];
D_temp = 0.02;
par.D = D_temp*ones(5,1); % generator rotational damping
par.M = [8; .1; 8; 8; .1];
D1 = par.D(1,1);
D2 = par.D(2,1);
D3 = par.D(3,1);
D6 = par.D(4,1);
D8 = par.D(5,1);
%D14 = par.D(2,1);
D = subs(D);
%%!!!50hz / 60hz
% generator inertias, as
syms M1 M2 M3 M6 M8 M14;
%TR250510 original
%M=(1/(pi*60))*[M1; M2; M3; M6; M8];
%original end
M=(1/(pi*50))*[M1; M3; M3; M6; M8]; %TR2502010 50Hz change
% (1/60*pi) * per unit "H"
% "H" values above
% (note: H is ratio of generator/turbine stored
% for gens @ buses % kinetic energy at rated speed,
% #’s 1,2,3,6,8 % to system wide "base" power
% (here chosen as 100 MVA). H is a standard
% spec in power system transient stability
% studies. It has units of seconds)
M1 = par.M(1,1);
M2 = par.M(2,1);
M3 = par.M(3,1);
M6 = par.M(4,1);
M8 = par.M(5,1);
%M14 = par.M(2,1);
M = subs(M);
HI = (deltime/2)*[diag(-(D)./M)-[[0;-(D(2:5))./M(2:5)] zeros(5,4)]...
diag(-ones(ngen+nslack,1)./M);...
[0 zeros(1,4); -ones(4,1) eye(ngen)] 0*eye(ngen+nslack)] ;
%%pfsolve calculates a first vbus estimation
% pfmiss
ibus=yb*vbus;
nmiss= vbus.*conj(ibus)+snet;
% pfmiss end
rmiss=[real(nmiss(genlist)); real(nmiss(loadlist));imag(nmiss(loadlist))];
error = max(abs(rmiss));
nsteps=5;
loop_continue = 1;
itcnt=0;
while (error > .00005 && loop_continue)
del=angle(vbus);
vmag=abs(vbus);
% pflowjac
ibus=yb*vbus;
dsdd=j*diag(conj(ibus).*vbus)-j*diag(vbus)*conj(yb)*diag(conj(vbus));
dsdv=diag(conj(ibus).*(vbus./abs(vbus)))+...
diag(vbus)*conj(yb)*diag(conj(vbus)./abs(vbus));
% end pflowjac
rjac = [ real(dsdd(genlist,genlist)) real(dsdd(genlist,loadlist)) ...
real(dsdv(genlist,loadlist)); real(dsdd(loadlist,genlist)) ...
real(dsdd(loadlist,loadlist)) real(dsdv(loadlist,loadlist)); ...
imag(dsdd(loadlist,genlist)) imag(dsdd(loadlist,loadlist))...
imag(dsdv(loadlist,loadlist))];
dx=- rjac\rmiss; %update
% update
ngen = length(genlist);%
nload = length(loadlist);%
nslack = length(slacklist);%
index1 = ngen;
index2 = index1 +1;
index3 = ngen + nload;
index4 = index3 + 1;
index5 = ngen + 2*nload;
del(genlist) = del(genlist) + dx(1:index1);
del(loadlist) = del(loadlist) + dx(index2:index3);
vmag(loadlist) = vmag(loadlist) + dx(index4:index5);
newvbus = vmag.*exp(j*del);
vbus = newvbus;
% END UPDATE
% pfmiss
ibus=yb*vbus;
nmiss= vbus.*conj(ibus)+snet;
% pfmiss end
rmiss=[real(nmiss(genlist)); real(nmiss(loadlist));imag(nmiss(loadlist))];
error = max(abs(rmiss));
itcnt= itcnt+1;
if itcnt >10
type convergefail.msg
loop_continue = 0;
end
end
xeq=[zeros(ngen+1,1); del([1;genlist])];
xs=0*xeq;% construct "natural" 10-dimensional
% state vector: generator freq deviations
% and generator phase angles. Note that
% freq deviation at equilibrium is forced to zero.
omega=0*ones(ngen+nslack,1); %omega vector initialisation
veq=vbus;
disp('The intial equilibrium power flow solution for the study network is:')
disp(' Bus# V p.u. Angle-degrees')
b_indices=(1:length(vbus));
disp([b_indices abs(vbus) (180/pi)*angle(vbus)])
%%initialisation values for embedded MATLAB function in simulink
ybinit=full(yb);
vbusinit=vbus;
omegainit=omega;
snetinit=snet;
xsinit=xs;
end

[Njuguna]

Its an initialisation code to embedded funtions.