nnodes = size(p, 2);      % Antal nodpunkter.

T = 2*pi;                 % Sluttid.

TimeSteps = 100;          % Antal tidssteg.

tider = 0:T/TimeSteps:T;  % Definierar tids-meshen.

xi = zeros(nnodes, TimeSteps+1);  % Initiera matris där lösningen lagras.
                                  % Varje kolonn innehåller lösningens
                                  % nodvärden vid en diskret tidpunkt.

				  
for i = 1:nnodes          % Loopa över alla nodpunkter.
  x = p(1, i);
  y = p(2, i);
  xi(i, 1) = u0_2D(x, y);
end

% Nu assemblerar vi matriser och vektorer som inte beror på tiden:

M = MinNya2DMassmatrisAssemblerare(p, t);  % Beräkna massmatrisen.

[A, Mc, R, rv] = MySecond2DPoissonAssembler(p, e, t);

% Nu startar vi tidsloopen:

for n = 1:TimeSteps  % Loopa över samtliga tidssteg.
  
  n                   % Skriver ut aktuellt tidssteg.
                      % Kan vara bra att veta hur
		      % långt programmet hunnit!
		      
  t_h = tider(n+1);   % Höger ändpunkt av tidsintervallet.
  t_v = tider(n);     % Vänster ändpunkt av tidsintervallet.
  dt = t_h - t_v;     % Tidsintervallets längd.
  
  b = Min2DLastAssemblerare(p, t, t_h);
  
  A_tot = A + Mc + R;
  b_tot = b + rv;
  
  xi_old = xi(:, n);    % Lösning vid tiden t = t_v.
  
  xi_new = (M + dt*A_tot)\(M*xi_old + dt*b_tot);  % Bakåt Euler!
  
  xi(:, n+1) = xi_new;  % Lösning vid tiden t = t_h.
  
end

pdesurf(p, t, xi(:, TimeSteps+1))
xlabel('x')
ylabel('y')
title('U(x,y,t), vid tiden t = T')
