Equazione del calore 2D

[Chebyshev1]

Buongiorno a tutti ^_^ Spero che qualcuno possa aiutarmi, perchè credo di avere delle lacune almeno sulla parte grafica e organizzativa di questo problema. Cercherò di spiegare passo passo come ho affrontato il problema.
Il problema da affrontare è:

$ (delT)/(delt) = k ((del^2T)/(delx^2)+(del^2T)/(dely^2)) $

Devo risolvere il problema su un dominio a L.

% Risoluzione dell'equazione del calore bidimensionale assegnate le
% condizioni al contorno di tipo lineare sui bordi NORD ed EST di un dominio
% a forma di L con passo spaziale uniforme e uguale lungo x e lungo x.
% E' stato utilizzato l'operatore del2 come approssimazione discreta
% dell'operatore differenziale del LAPLACIANO.
%
%Si ricordi che:
% l = del2(phi) = (1/4)*∆ = (1/4)*[(d²phi/dx²)+(d²phi/dy²)]
% e che quindi
% l(ij) = (1/4)*[(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1))]-phi(i,j)
% sui nodi interni, mentre sui nodi di bordo è applicata un'interpolazione
% cubica.
%
% Il reticolo di nodi utilizzato è del tipo:
%
%
%
% --------------------------------------------> x
% |
% | x -- x -- x -- x -- x -- x -- x -- x
% |
% | x -- o -- o -- o -- o -- o -- o -- x
% |
% | x -- o -- o -- o -- o -- o -- o -- x
% |
% | x -- o -- o -- o -- o -- o -- o -- x
% |
% | x -- x -- x -- x -- o -- o -- o -- x
% |
% | x -- x -- x -- x -- o -- o -- o -- x
% |
% | x -- x -- x -- x -- o -- o -- o -- x
% y |
% ˅
%
%% Dati di input
k = 1.0; % Coefficiente di diffusione
Bet = 0.245; % Parametro diffusivo
lx = 1.0; % Lunghezza del dominio lungo la direzione x
ly = 1.0; % Lunghezza del dominio lungo la direzione y
nx = 8; % numero di nodi lungo la direzione x
ny = 8; % numero di nodi lungo la direzione y
N = 20; % numero di nodi temporali
%
%% Dati ricavabili da quelli di input
Dx = lx/(nx-1); % passo spaziale lungo x
Dy = ly/(ny-1); % passo spaziale lungo y
x = 0:Dx:lx; % Mesh spaziale dei nodi in direzione x
y = 0:Dy:ly; % Mesh spaziale dei nodi in direzione y
l = lx; % = ly
n = nx; % = ny
h = Dx; % = Dy
Dt = Bet*(h^2)/k; % Passo temporale
t = 0:Dt:N*Dt; % Mesh temporale
%
%% Costruzione della griglia
G=numgrid('L',n);
%
%% Condizioni al contorno alla Dirichlet
phi0=0*G;
phi0(1,:)=linspace(0,1,n); %condizione lineare sul contorno NORD
phi0(:,end)=linspace(1,0,n); %condizione lineare sul contorno EST
%
%% Costruzione del Laplaciano secondo
LB=4*del2(phi0); %Laplaciano secondo
%Ricordiamo che le nostre incognite sono i valori nodali; nella risoluzione
%dell'equazione differenziale calcoleremo il del2 ovunque, dobbiamo tener
%presente, però, che quello calcolato al contorno è stato ottenuto con
%un'interpolazione cubica e quei valori non sono di nostro interesse,
%perchè sono già noti.
L=LB(G>0); %è il vettore contenente le componenti del termine noto del
%Laplaciano secondo a meno dei parametri di griglia
%
%% Elimino dalla condizione iniziale i valori di bordo
phiold=phi0(G>0);
%
%% Ragiono sulle soluzioni valutate come vettori colonna
neq=max(G(:));
phinew=zeros(neq,N+1);
phi(:,:,1)=phi0;
for it=2:N+1
phinew(:,it)=phiold+k*Dt*L;
%quello che ho ottenuto è un vettore colonna che devo riordinare all'interno della mia griglia
K=logical(G);
s=zeros(n);
s(K)=phinew(:,it); %in questo modo sono ordinati nella griglia secondo il nostro schema
phiold=phinew(:,it);
phi(:,:,it)=phi0+s;
end

Con questo che ho fatto ho creato N+1 matrici, che rappresentano quelle al variare del tempo della mia soluzione sui nodi del problema.
Credete sia corretto? Ma soprattutto: come faccio a rappresentare tutto ciò graficamente?
Grazie mille a chi mi dedicherà del tempo

[Raptorista1]

Non ho letto il codice a fondo, ma mi sembra di capire che usi differenze finite in spazio e tempo [eulero esplicito].

Per la rappresentazione grafica, probabilmente la cosa più semplice è usare la funzione contour di matlab/octave, che fa un grafico facilmente interpretabile di una matrice di valori.