  
% Equazione del calore
c=0.5;
% discretizzazione nello spazio
sigma=0; % consideriamo solo il caso "sigma=0"
N=18;     % numero di nodi della griglia
% qui serve anche il passo della griglia (viene usato per verificare la condizione di stabilità di Eulero in avanti)
h=1/(N-1);     % passo della griglia

% costruisco la matrice "A" (approssimazione del laplaciano):
A = mat_ellipt_2D(N,sigma);
%A=-A;
u_vera=ones(size(A,2),1); %voglio che la soluz esatta sia questa
b=A*u_vera; 

x=[h:h:1-h]';
y=[h:h:1-h]';
[xx,yy]=meshgrid(x,y);

esperimento=0;  % 0=riscaldamento,  1=raffreddamento

% impongo le condizioni iniziali
if esperimento==0
  u0=zeros(size(A,2),1);       % esempio di "riscaldamento di un solido"
else
  u0=zeros(size(A,2),1)+2;   % esempio di "raffreddamento di un solido"
end

%
% discretizzazione nel tempo:  metodo di Eulero esplicito
%
figure(1); clf, 
num=1;
t0=0;
T=0.2;
dt = 0.25 * (h^2 / (2*c));  % scelgo il passo di discretizzazione temporale
tt=[t0:dt:T];
u=u0;
for t=tt(2:end-1)
%    % costruisco "b":
%    b=zeros(dim,1);
%    k=1;
%    for j=1:(sqrt(n)-2)
%       for i=1:(sqrt(n)-2)
%           %b(k)=sin(x(i)*y(j))-cos(t);     % ad esempio con questa funzione
%           %b(k)=0; % esempio di "raffreddamento di un solido"
%           b(k)=10; % esempio di "riscaldamento di un solido"
%           k=k+1;
%       end
%    end
%    % metodo di Eulero esplicito
    v=u(:,end);
    v = v + dt*(-A*v+b);
    u=[u v];
    % Grafico della soluzione in "t" di cui ho calcolato il valore sui nodi della griglia
    sol_t=reshape(u(:,end),N-2,N-2);
    mesh(xx,yy,sol_t); 
    if esperimento==0  axis([0 1 0 1 0 1]); else  axis([0 1 0 1 0 2]); end; 
    title(['t = ' num2str(t)]), drawnow, tic, 
    while toc<0.05, end;
end

