Matlab - Tengo un problema en la programacion del metodo de los paneles en una placa plana

Hola a todos les comento que tengo que hacer un programa, sobre fuerzas de sustentacion en perfiles delgado de envergadura finita. Utilizo el Metodo de los paneles, que basicamente consiste en utilizar la ley de Biot Sabart donde en vez del campo magnetico inducido es una velocidad. Por condicion de impermeabilidad esta velocidad en el punto central del panel debe ser cero.
Aqui los dejo el programa que hice pero no esta funcionando por la matriz que llamo Wm no queda diagonal dominante.
Se puede ver que el programa no devuelve resultados admisibles cuando vemos los graficos de circulacion segun envergadura o cuerda porque la circulacion en los extremos tiene que tender a cero.

Espero que alguien me pueda ayudar porque ya no se que mas hacer, hice el programa de la manera que sea mas facil para la lectura, cualquier duda estoy a su disposicion.

Saludos


% Vela de envergadura finita
% Metodo de los paneles
% El plano cuerda-envergadura se encuentra en el plano x-z
%
% Ingrese las dimensiones de la vela
puj=input('Teclee la longitud del pujamen (en metros):');
gra=input('Teclee la longitud del gratil (en metros):');
tope=input('Teclee la longitud del tope (en metros):');
%
% Ingrese la cantidad de divisiones para generar el mallado
bas=input('Teclee cantidad de nodos en la base:');
altu=input('Teclee cantidad de nodos en el gratil:');
%
% Ingresamos los datos de la corriente externa para saber las condiciones
% de contorno
U=input('Teclee la velocidad de corriente libre en m/s:');
ro=input('Teclee la densidad del aire en kg/m^3:');
al=input('Teclee el angulo de ataque en grados:');
% la matriz zz es una matriz de prueba que genera con los datos que se
% ingresa las coordenadas z de los puntos
% La matriz y tiene las coordenadas y de cada nodo del mallado
y=[];
for k=1:bas
for h=1:altu
y(k,h)=(gra/(altu-1))*h-(gra/(altu-1));
end
end
y;
% La matriz xx y yy son matrices para generar la matriz x por medio de
% interpolacion lineal
xx=[0:puj/(bas-1):puj;0:tope/(bas-1):tope]';
yy=[0:gra/(altu-1):gra];
x=[];
for k=1:bas
for h=1:altu
if k==1 && h==1 && k==bas && h==altu
x(k,h)=xx(k,h);
else
x(k,h)=xx(k,1)+((xx(k,2)-xx(k,1))/(gra))*yy(h);
end
end
end
x;
figure(1),mesh(x,y,zeros(bas,altu)),title('Grafico de los paneles')...
,xlabel('x[m]'),ylabel('y[m]'),axis equal
desp=[];
for k=1:bas
for h=1:altu
if k~=bas
desp(k,h)=0.25*(x(k+1,h)-x(k,h));
else
desp(k,h)=desp(k-1,h);
end
end
end
% Para las posiciones de los nodos de los paneles son la matriz z y la
% matriz X=x+desp, y Y es para las pocisiones en y
X=x+desp;
Y=y;
% Uso el comando global para que las pueda ver en las subrutinas tambien
X;
Y;
figure(2),surf(X,Y,zeros(bas,altu)),title('Grafico de los anillos vorticosos')...
,xlabel('x[m]'),ylabel('y[m]'),axis equal
% Ahora hay que calcular los centros de los paneles utilizandon las
% matrices X y Y
% Las coordenadas de los centros se almacenan en Xm y Ym (puntos de
% control)
for k=1:(bas-1)
for h=1:(altu-1)
Xm(k,h)=(X(k+1,h)-X(k,h))/2+X(k,h);
Ym(k,h)=(Y(k,h+1)-Y(k,h))/2+Y(k,h);
end
end
Xm;
Ym;
disp('Se termino de calcular los puntos de control')
%--------------------------------------------------------------------------
% RECORDAR ESTE MODULO SIRVE UNICAMENTE SI LA VELA ES PLANA CONSIDERO UNA
% UNICA NORMAL
% En caso de haber curvatura hay que generalizar a la normal y va a tener
% mas componentes
% N tiene componentes en Z unicamente
% La matriz A contiene los coeficientes del sistema a resolver Wm,n *
% (gamma)=Wm
xm=0;
ym=0;
cont=0;
cont2=0;
Wm=zeros((bas-1)*(altu-1),(bas-1)*(altu-1));
for q=1:(bas-1)
for p=1:(altu-1)
xm=Xm(q,p);
ym=Ym(q,p);
cont2=cont2+1;
cont=0;
w=zeros(1,(bas-1)*(altu-1));
for k=1:(bas-1)
for h=1:(altu-1)
cont=cont+1;
if k==(bas-1)
% Por aca anda el quilombo
w(cont)=((1/((xm-X(k,h))*(ym-Y(k,h+1))-(xm-X(k,h+1))*(ym-Y(k,h))))*(((X(k,h+1)-X(k,h))*(xm-X(k,h))+((Y(k,h+1)-Y(k,h))*(ym-Y(k,h))))/(sqrt((xm-X(k,h))^2+(ym-Y(k,h))^2))-...
((X(k,h+1)-X(k,h))*(xm-X(k,h+1))+(Y(k,h+1)-Y(k,h))*(ym-Y(k,h+1)))/(sqrt((xm-X(k,h+1))^2+(ym-Y(k,h+1))^2)))+...
(1/(Y(k,h)-ym))*(1+(xm-X(k,h))/(sqrt((xm-X(k,h))^2+(xm-Y(k,h))^2)))-...
(1/(Y(k,h+1)-ym))*(1+(xm-X(k,h+1))/(sqrt((xm-X(k,h+1))^2+(xm-Y(k,h+1))^2))));
else
r0a=[X(k,h+1)-X(k,h) Y(k,h+1)-Y(k,h) 0];
r1a=[xm-X(k,h) ym-Y(k,h) 0];
r2a=[xm-X(k,h+1) ym-Y(k,h+1) 0];
r0b=[X(k+1,h+1)-X(k,h+1) Y(k+1,h+1)-Y(k,h+1) 0];
r1b=[xm-X(k,h+1) ym-Y(k,h+1) 0];
r2b=[xm-X(k+1,h+1) ym-Y(k+1,h+1) 0];
r0c=[X(k+1,h)-X(k+1,h+1) Y(k+1,h)-Y(k+1,h+1) 0];
r1c=[xm-X(k+1,h+1) ym-Y(k+1,h+1) 0];
r2c=[xm-X(k+1,h) ym-Y(k+1,h) 0];
r0d=[X(k,h)-X(k+1,h) Y(k,h)-Y(k+1,h) 0];
r1d=[xm-X(k+1,h) ym-Y(k+1,h) 0];
r2d=[xm-X(k,h) ym-Y(k,h) 0];
w(cont)=(dot(r0a/(norm(cross(r1a,r2a))),(r1a/norm(r1a)-r2a/norm(r2a)))*(dot(cross(r1a,r2a),[0 0 1]))/norm(cross(r1a,r2a))+...
dot(r0b/(norm(cross(r1b,r2b))),(r1b/norm(r1b)-r2b/norm(r2b)))*(dot(cross(r1b,r2b),[0 0 1]))/norm(cross(r1b,r2b))+...
dot(r0c/(norm(cross(r1c,r2c))),(r1c/norm(r1c)-r2c/norm(r2c)))*(dot(cross(r1b,r2b),[0 0 1]))/norm(cross(r1c,r2c))+...
dot(r0d/(norm(cross(r1d,r2d))),(r1d/norm(r1d)-r2d/norm(r2d))))*(dot(cross(r1b,r2b),[0 0 1]))/norm(cross(r1d,r2d));
end
end
end
Wm(cont2,:)=w(:);
end
end
Wm;
disp('Se calculo la matriz de coeficientes')
% Caluculamos la inversa de Wm para resolver el sistema de ecuaciones
% lineales
Wminv=inv(Wm);
disp('Se calculo la inversa de la matriz de coeficientes')
%--------------------------------------------------------------------------
% w es el vector de terminos independientes
w=zeros((altu-1)*(bas-1),1);
alpha=al*pi/180;
for k=1:((altu-1)*(bas-1))
w(k,1)=-4*pi*alpha*U;
end
w;
disp('Se calculo el vector independiente')
Gamma=Wminv*w;
disp('Se resolvio el sistema')
% Creamos una matriz gm similar a Xm e Ym donde se almacenen las circulaciones
% de cada panel
cont3=0;
gm=zeros((bas-1),(altu-1));
for k=1:(bas-1)
for h=1:(altu-1)
cont3=cont3+1;
gm(k,h)=Gamma(cont3);
end
end
% Matriz para calcular las circulaciones segun cuerda Calculamos las
% guerzas de sustentacion individuales
%
Gm=zeros(bas-1,altu-1);
l=zeros((bas-1),(altu-1));
for k=1:(bas-1)
for h=1:(altu-1)
if k==1
l(k,h)=-ro*U*gm(k,h);
Gm(k,h)=gm(k,h);
else
l(k,h)=-ro*U*(-gm(k-1,h)+gm(k,h));
Gm(k,h)=(-gm(k-1,h)+gm(k,h));
end
end
end
l;
L=0;
G=0;
% Calculamos la suma de las fuerzas de sustentacion de cada elemento
for k=1:(bas-1)
for h=1:(altu-1)
L=L+l(k,h);
G=G+gm(k,h);
end
end
L
A=(tope+puj)*gra/2;
Cl=L/(0.5*ro*(U^2)*A)
% Grafico de variacion de circulacion segun envergadura
for k=1:(bas-1)
figure(3),plot(Ym(k,:)/gra,Gm(k,:)),title('Grafico de variacion de Circulacion segun envergadura')...
,xlabel('y/gra'),ylabel('\Gamma*'),grid
hold on
end
hold off
% Grafico de circulacion segun envergadura
for k=1:(bas-1)
figure(4),plot(Ym(k,:)/gra,gm(k,:)),title('Grafico de Circulacion segun envergadura')...
,xlabel('y/gra'),ylabel('\Gamma*'),grid
hold on
end
hold off
% Grafico de circulacion segun cuerda
for h=1:(altu-1)
figure(5),plot(Xm(:,h)/puj,gm(:,h)),title('Grafico de Circulacion segun cuerda')...
,xlabel('x/pujamen'),ylabel('\Gamma*'),grid on
hold on
end
hold off
% Grafico de variacion circulacion segun cuerda
for h=1:(altu-1)
figure(6),plot(Xm(:,h)/puj,Gm(:,h)),title('Grafico de variacion Circulacion segun cuerda')...
,xlabel('x/pujamen'),ylabel('\Gamma*'),grid on
hold on
end
hold off

Hola Martín, ¿conseguiste resolver tu problema? Es que quiero tratar de hacer lo mismo, es decir, hacer el método de paneles en matlab y me gustaría ponerme en contacto contigo.