Matlab - Problema con Funcion de tensiones en barra eliptica/rectangular

Hola!,

un gusto a todos... estoy teniendo inconvenientes para poder realizar la funcion en Matlab de forma correcta.

Nose donde estoy teniendo la complicacion.

Podran ayudarme?

%************************************************************************%
% G:modulo de elasticidad transversal%
% A:angulo de torsion especifico por unidad de longitud%
% a:eje horizontal de la seccion eliptica%
% b:eje vertical de la seccion eliptica %
% tol:tolerancia%
% max:maxima cantidad de veces que es recorrida la malla%
% xi: coordenada en el eje x del punto a calcular%
% yi: coordenada en el eje y del punto a calcular%
%************************************************************************%
function U=torsion_barra_eliptica_(G,A,a,b,h,tol,max1,xi,yi)
n=fix((a/h)+0.0000000001)+1;
m=fix((b/h)+0.0000000001)+1;
c=a./2;
d=b./2;
Vi=0;
U=Vi*ones(n,m);
w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));
err=1;
cnt=0;
while((err>tol)&(cnt<=max1))
err=0;
for j=2:m-1
for i=2:n-1
relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4.*U(i,j)-(h.^2).*(-2.*G.*A))/4;
U(i,j)=U(i,j)+relx;
if (err<=abs(relx))
err=abs(relx);
end
x=h.*(i-(n+1)./2);
y=h.*(j-(m+1)./2);
if ((x.^2)./(c.^2)+(y.^2)./(d.^2)>1);
U(i,j)=0;
end
end
end
cnt=cnt+1;
end
U=flipud(U');
W=U;
x=(-a./2):h:(a./2);
y=(-b./2):h:(b./2);
f=W(fix(((y+(b./2))./h)+1+0.0000000001),fix(((x+(a./2))./h)+1+0.0000000001));
[fx fy]=gradient(f,h);
[x y]=meshgrid(x,y);
subplot(2,2,1)
mesh(x,y,f)
xlabel('X')
ylabel('Y')
zlabel('Fi')
title('FUNCION DE TENSIONES')
subplot(2,2,2)
quiver(x,y,fx,fy)
title('CAMPO DE GRADIENTE')
Fx=fy;
Fy=-fx;
subplot(2,2,3)
quiver(x,y,Fx,Fy)
title('CAMPO DE TENSIONES TANGENCIALES')
subplot(2,2,4)
contour(f,20)
title('LINEAS DE NIVEL DE LA FUNCION DE TENSIONES')
M=sqrt((fx(fix((yi+(b./2))./h+1),fix((xi+(a./2))./h+1))).^2+(fy(fix((yi+(b./2))./h+1),fix((xi+(a./2))./h+1))).^2);
disp('La tension en el punto (xi,yi) es:')
disp(M)