Script with Klikenberg effect

can anyone help me, why my script does not work correctly. the MatLab can not proses the function "Klinkenberg"

and here is the script,

Pm=3800; % Initial Reservoir Pressure

% length of block

dx = 22;

dy = 12;

dz = 2;

dt = 10; % Time Interval

T = 1000; % Total days

n = 5; % no. of matrix

i=1; j=1; k=1;

Ax=dy*dz;

Ay=dx*dz;

Az=dx*dy;

B(n,n,n) = 0;

S(n,n,n) = 0;

W(n,n,n) = 0;

E(n,n,n) = 0;

N(n,n,n) = 0;

A(n,n,n) = 0;

X(n,n,n) = 0;

m=T/dt;

P(n,n,n,m) = 0;

dP(n,n,n) = 0;

Km = klinkenberg(Pm);

Tgsc= trans(Pm);

Sgm=1;

B(:,:,:)= (Km * Az*Tgsc)/dz;

S(:,:,:)= (Km * Ay*Tgsc)/dy;

W(:,:,:)= (Km * Ax*Tgsc)/dx;

E(:,:,:)= (Km * Ax*Tgsc)/dx;

N(:,:,:)= (Km * Ay*Tgsc)/dy;

A(:,:,:)= (Km * Az*Tgsc)/dz;

X(:,:,:)= ((1)*(dx*dy*dz)/dt)*((Sgm*mat_por(Pm)*cmprs(Pm) *rho_sc()/(5.61458*Bg_factor(Pm)))+( mat_por(Pm)*cmprs(Pm)*(ad_vol(Pm)+de_vol(Pm))));

Q = -1*X*Pm;

C = (E+W+N+S+A+B-X);

tun(125,125)=0;

lun(125,1)=0;

for y=1:int16(fix((T/dt)))

for x=1:125

[i,j,k]=index(x);

lun(x,1)=Q(i,j,k);

temp=(25*(i-1))+(5*(j-1))+k;

tun(x,temp)= C(i,j,k);

i1 = i-1;

if i1>=1

temp=(25*(i1-1))+(5*(j-1))+k;

tun(x,temp)=W(i,j,k);

end

j1 = j-1;

if j1>=1 temp=(25*(i-1))+(5*(j1-1))+k;

tun(x,temp)=S(i,j,k);

end

k1 = k-1;

if k1>=1

temp=(25*(i-1))+(5*(j-1))+k1;

tun(x,temp)=B(i,j,k);

end

i1=i+1;

if i1<=5

temp=(25*(i1-1))+(5*(j-1))+k;

tun(x,temp)=E(i,j,k);

end

j1=j+1;

if j1<=5

temp=(25*(i-1))+(5*(j1-1))+k;

tun(x,temp)=N(i,j,k);

end

for i=1:n

for j=1:n

for k=1:n

if y==1

P(i,j,k,y)= Pm -dP(i,j,k);

else P(i,j,k,y)= P(i,j,k,y-1)-dP(i,j,k);

end

end

end

end

for i=1:n

for j=1:n

for k=1:n

Km=klinkenberg( P(i,j,k,y) );

Tgsc = trans( P(i,j,k,y) );

B(i,j,k)= (Km * Az*Tgsc)/dz;

S(i,j,k)= (Km * Ay*Tgsc)/dy;

W(i,j,k)= (Km * Ax*Tgsc)/dx;

E(i,j,k)= (Km * Ax*Tgsc)/dx;

N(i,j,k)= (Km * Ay*Tgsc)/dy;

A(i,j,k)= (Km * Az*Tgsc)/dz;

X(i,j,k)= ((1)*(dx*dy*dz)/dt)*((Sgm*mat_por(P(i,j,k,y))*cmprs(P(i,j,k,y))*rho_sc()/(5.61458 *Bg_factor(P(i,j,k,y))))+(mat_por(P(i,j,k,y))*cmprs(P(i,j,k,y))*(ad_vol(P(i,j,k,y))+ de_vol(P(i,j,k,y)))));

Q(i,j,k) = -1*X(i,j,k)*P(i,j,k,y); C(i,j,k) = (E(i,j,k)+W(i,j,k)+N(i,j,k)+S(i,j,k)+A(i,j,k)+B(i,j,k)-X(i,j,k));

end

end

end

end

P1(n,n,n,m)=0;

for y=1:m

for i=1:n

for j=1:n

for k=1:n

P1(i,j,k,y)=P(k,i,j,y);

end

end

end

end