index error problem how can fix it?

This is the code that what i put in but there is a some error the line at abc2(a_tic)=Vexact(xind,yind,zind)

error says that Index at position 3 exceeds array boundaries (must not exceed 11).

but i have no idea to fix it

how can i do it?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FINITE DIFFERENCE SOLUTION OF POISSON'S EQUATION:

% Vxx+Vyy+Vzz=G

% USING THE METHOD OF SUCCESSIVE OVER-RELAXATION

%

% NX : NO.OF INTERVALS ALONG X-AXIS

% NY : NO.OF INTERVALS ALONG Y-AXIS

% NZ : NO.OF INTERVALS ALONG Z-AXIS

%A X B : DIMENSION OF THE SOLUTION REGION

% V(I,J,K) : POTENTIALAT GRID POINT(X,Y Z)=H*(I,J,K)

% WHERE I=0,1,......,NX,J=0,1,......,NY K=0,1,......,NZ

% H :MESH SIZE

%***************************************************************

%SPECIFY BOUNDARY VALUES TO ZEROS OFR TO AVERAGE OF FIXED VALUES

A=1; B=1; C=1;

V1=0; V2=10; V3=20; V4=-10; V5=20; V6=40;

NX=20; % 4 12 20

NY=NX;

NZ=NX;

H=A/NX;

%SET INITIAL GUESS EQUAL TO ZEROS OR TO AVERAGE OF FIXED VALUES

for I=1:NX-1

for J=1:NY-1

for K=1:NZ-1

V(I+1,J+1,K+1)=(V1+V2+V3+V4+V5+V6)/6.0;

end

end

end

%SET POTENTIALS AT FIXED NODES

for K=1:NZ-1

for J=1:NY-1

V(1,J+1,K+1) = V1;

V(NX+1,J+1,K+1)=V4;

end

end

for I=1:NX-1

for K=1:NZ-1

V(I+1,1,K+1)=V2;

V(I+1,NY+1,K+1)= V5;

end

end

for I=1:NX-1

for J=1:NY-1

V(I+1,J+1,1)=V3;

V(I+1,J+1,NZ+1)=V6;

end

end

V(1,1,1)=(V1+V2+V3)/3.0;

V(NX+1,1)=(V1+V4)/2.0;

V(1,NY+1,1)=(V2+V5)/2.0;

V(1,1,NZ+1)=(V3+V6)/2.0;

V(NX+1,NY+1,NZ+1)=V4+V5+V6/3.0;

%FIND THE OPTIMUM OVER-RELAXATION FACTOR

T=cos(pi/NX)+cos(pi/NY)+cos(pi/NZ);

W=(8-sqrt(64-16*T^2))/(T^2);

disp(['SOR Factor Omega=',num2str(W)])

W6=W/6;

%INERATION BEGINS

NCOUNT = 0;

loop= 1;

while loop == 1;

RMIN =0;

for I=1:NX-1

X=H*I;

for J=1:NY-1

Y=H*J;

for K=1:NZ-1

Z=H*K;

G=-36.0*pi*X*(Y-1.0);

R=W6*(V(I+2,J+1,K+1)+V(I,J+1,K+1)+V(I+1,J+2,K+1)+V(I+1,J,K+1)+V(I+1,J+1,K)-6.0*V(I+1,J+1,K+1)-G*H*H*H);

RMIN=RMIN+abs(R);

V(I+1,J+1,K+1)=V(I+1,J+1,K+1)+R;

end

end

end

RMIN=RMIN/(NX*NY+NZ);

if(RMIN>=0.0001)

NCOUNT=NCOUNT +1;

if(NCOUNT>100)

loop=0;

disp('SOLUTION DOES NOT CONVERGE IN 100 INTERATIONS')

end

else %Then RMIN is less than.0001 and then solution has converged

loop=0;

disp(['Solution Converges in',num2str(NCOUNT),'iterations'])

end

end

Vnum=V;

%Grab original points a through i

abc=zeros(1,9);

a_tic=1;

vec=[0:H:1];

for ii=.25:.25:.75

for jj=.25 :.25 :.75

for kk=.25:.25:.75

xind=find(vec==ii);

yind=find(vec==jj);

zind=find(vec==kk);

%disp([xind,yind,zind])

abc(a_tic)=Vnum(xind,yind,zind);

a_tic=a_tic+1;

end

end

end

%OUTPUT THE FINITIE DIFFERENCE APPROX.RESULTS

%----------------------------------------------------------------

%CALCULATE THE EXACT SOLUTION

%

%POISSON'S EQUATION WITH HOMOGENEOUS BOUNDARY CONDITIONS

%SOLVED BY SERIES EXPANSION

for I=1:NX-1

X=H*I;

for J=1:NY-1

Y=H*J;

for K=1:NZ-1

Z=H*K;

end

end

end

SUM=0;

for M=1:10 % TAKE ONLY 10 TERMS OF THE SERIES

FM=M;

for N=1:10

FN=N;

for W=1:10

FW=W;

FACTOR1=(FM*pi/A)^2 +(FN*pi/B)^2+(FW*pi/C)^2;

FACTOR2=((-1)^(M+N+W))*(-144)*A*B*C^2/(pi^2*FM*FN*FW);

FACTOR3=1+1/W*pi;

FACTOR=FACTOR2*FACTOR3/FACTOR1;

SUM=SUM+FACTOR*sin(FM*pi*X/A)*sin(FN*pi*Y/B)*sin(FW*pi*Z/C);

end

end

end

VH=SUM;

%LAPLACE'S EQUATION WITH INHOMOGENEOUS BOUNDARY CONDITIONS

%SOLVED USING THE METHOD OF SEPARATION OF VARIABLES

C1=16*V1/pi^2;

C2=16*V2/pi^2;

C3=16*V3/pi^2;

C4=16*V4/pi^2;

C5=16*V5/pi^2;

C6=16*V6/pi^2;

SUM=0;

for K=1:10 %TAKE ONLY 10 TERMS OF THE SERIES

N=2*K-1;

for M=2*K-1;

for W=2*K-1;

AN=N;,AM=M;,AW=W;

A1=sin(AN*pi*Y/B);

A2=sin(AW*pi*Z/C);

A3=sinh(AN*AW*pi*X/B*C);

A4=AN*AW*sinh(AN*AW*A*pi/B*C);

TERM1=C1*A1*A2*A3/A4;

B1=sin(AM*pi*X/A);

B2=sin(AW*pi*Z/C);

B3=sinh(AM*AW*pi*Y/A*C);

B4=AM*AW*sinh(AM*AW*B*pi/A*C);

TERM2=C2*B1*B2*B3/B4;

D1=sin(AM*pi*X/A);

D2=sin(AN*pi*Y/B);

D3=sinh(AM*AN*pi*Z/A*B);

D4=AM*AN*sinh(AM*AN*C*pi/A*B);

TERM3=C3*D1*D2*D3/D4;

E1=sin(AN*pi*Y/B);

E2=sin(AW*pi*Z/C);

E3=sinh(AN*AW*pi*(A-X)/B*C);

E4=AN*AW*sinh(AN*AW*A*pi/B*C);

TERM4=C4*E1*E2*E3/E4;

F1=sin(AM*pi*X/A);

F2=sin(AW*pi*Z/C);

F3=sinh(AM*AW*pi*(B-Y)/A*C);

F4=AM*AW*sinh(AM*AW*B*pi/A*C);

TERM5=C5*F1*F2*F3/F4;

G1=sin(AM*pi*X/A);

G2=sin(AN*pi*Y/B);

G3=sinh(AM*AN*pi*(C-Z)/A*B);

G4=AM*AN*sinh(AM*AN*C*pi/A*B);

TERM6=C6*G1*G2*G3/G4 ;

TERM=TERM1+TERM2+TERM3+TERM4+TERM5+TERM6;

SUM=SUM+TERM;

end

end

end

VI=SUM;

Vexact(I+1,J+1,K+1)=VH+VI;

%Grab original points a through i

abc2=zeros(1,9);

a_tic=1;

vec=[0:H:1];

for ii=.25:.25:.75

for jj=.25:.25:.75

for kk=.25:.25:.75

xfind=find(vec==ii);

yfind=find(vec==jj);

zfind=find(vec==kk);

%disp([xind,yind,zind])

abc2(a_tic)=Vexact(xind,yind,zind) <------------------------------------------here is the error region

a_tic=a_tic+1;

end

end

end

figure(1);

imagesc(flipud(Vnum')),

colorbar

xlabel('x'), ylabel('y'),zlabel('z')

title('Example 3.4 : Poisson PDE')

format short g

[G A]

May be you have to rename xfind, yfind and zfind in the loop with error as xind,yind and zind? Or use

abc2(a_tic)=Vexact(xfind,yfind,zfind) 

[ljaseon]

Thanks for help i got it