[Matlab]Risoluzione con Fourier eq.derivate parziali pura diffusione

[entrochim]

Buonasera a tutti,
ho un problema con Matlab nella stesura del codice per la risoluzione di una equazione a derivate parziali parabolica in un dominio bidimensionale. Il tutto da risolvere con Fourier.
Ut=A(Uxx+Uyy)
In pratica si tratta di un problema di pura diffusione di calore all'interno di un'aletta di uno scambiatore di calore con analisi sia del transitorio che dello stato stazionario.
Link all'immagine per un'idea della geometria: http://i.imgur.com/xljuvZf.png
L'equazione è stata adimensionalizzata e le condizioni al contorno rese omogenee. In particolare a sinistra abbiamo condizioni di Dirichlet, a destra Neumann, sopra e sotto miste (Robin). Le condizioni miste complicano di molto la risoluzione perchè si ottiene un'equazione del tipo cotan(x)=(Bi^2-x^2)/(Bi*x)
Allego la risoluzione: http://i.imgur.com/PQhjpzA.png

Lo scheletro di codice che ho fatto è il seguente

W = 10;             %larghezza
H = 10;             %altezza
T = 20;         %T(x,y,0) ambiente al tempo zero
tfinal=50;  % tempo finale
 
Npt=50;           %Nodi
x=linspace(0,W,Npts);
y=linspace(0,H,Npts); %Nodi in questo caso è quadrato ma vabbè ci vuol poco a renderlo rettangolare
[X,Y]=meshgrid(x,y);


%Parte dove valuto gli autovalori sull'asse delle y con le condizioni miste omogenee di green. 'q' è l'autovalore, c il numero di Biot

for e=1:Npt
qfun=@(q,c)cot(q)-(q^2-c^2)/(2*c*q);
c= 2;
fun=@(q)qfun(q,c);
r = fzero(fun,e);
end

%A questo punto devo trovare i valori dei miei coefficienti anm. Ci va messa quella bella lunga relazione con integrale doppio sopra e sotto della seconda immagine.

for jj=1:Npt
amn = La relazione di prima con integrale dxdy diviso int della roba al quadrato
   a(jj) = amn;
end

%Trovati tutti i coefficienti posso passare al ciclo finale dove valuto tutte le temperature nei vari nodi. Imposto la serie di Fourier.

U=zeros(Npts,Npts,t);
for t=1:tfinal
	for jj=1:Npt       
  	Un=a(jj)(exp(-(r+ sin((jj-1/2)*pi*x/W)))*sin((jj-1/2)*pi*x/W).*((cos(r))+ c*sin(r)/ r);
 	 U=U+Un;
	end
end

Parte dei grafici da vedere

L'impostazione è corretta? Trovare prima i coefficienti della serie di fourier e poi risolvere la serie stessa? I cicli for vanno separati o possono essere accorpati?
Il ciclo for per valutare l'autovalore delle condizioni miste ha un senso?