L = 1;
h = 0.1;
x = [0:h:L];
N = length(x);
k=1;
b=0; % coefficiente di convezione
beta = -20;  % coefficiente di reazione
f = zeros(1,N);
% definisco la funzione rilevamento del bordo:
g = zeros(N,1);
g(1) = 0.9;
g(N) = 0.9;
% calcolo le matrici del modello discretizzato nello spazio:
[M,K,F] = calc_discr_Galerkin_1D(N,diff(x),k,b,f);
% applico le condizioni al contorno di Dirichlet:
F(2) = F(2) - g(1)*(K(2,1)+beta*M(2,1));
F(N-1) = F(N-1) - g(N)*(K(N-1,N)+beta*M(2,1));
B=zeros(N,N-2);
for i=2:N-1  B(i,i-1)=1; end;
M = M(2:N-1,2:N-1);
K = K(2:N-1,2:N-1);
F = F(2:N-1);
%
% definisco la condizione iniziale:
u0 = zeros(N,1);
u0(1) = 0.9;
u0(N) = 0.9;
%
Tf = 0.03;
dt = 0.001;
% Analisi della matrice di iterazione:
rho_eig = sort(eig(full(inv(M)*K)))
sort(dt*(beta+rho_eig))
autov_mat_iteraz = sort(eig((1-dt*beta)*eye(size(M))-dt*inv(M)*K))
%
u_nd = u0(2:N-1);
for i=1:floor(Tf/dt)
  u_nd = M \ (((1-dt*beta)*M-dt*K)*u_nd + dt*F);
  u = B*u_nd + g;
  figure(1); plot(x,u,'b.-'); axis([0 L -0.1 1]); title(['t = ' num2str(i*dt)]), drawnow, pause(0.3);
end