function U = MyHeatEqSolver(p, t, e, InitialData, EqData, BdryData, FinalTime)

% Solves the heat equation cdu/dt-D(aDu)=f in Omega, with heat capacity c and
% heat conductivity a, and with bdry conditions dudn = k(g-u) and initial
% condition u = u0.
% The data a, c and f are specified in EqData, the bdry data k and g are
% specified in BdryData, and u0 is specified in InitialData (see the example
% code below). The initial condition corresponds to time = 0 and the solution
% is computed up to FinalTime. The geometry and spatial elements are specified
% in p, t and e (as in MyPoissonSolver). The solution at time j*k is stored in
% U(:,j+1).

time = 0;

u = feval(InitialData, p);
if size(u,1)==1 u=u'; end

U = u;

k = 0.01;

h = waitbar(0, 'Just a minute .. ');

while time < FinalTime

  time = time + k/2;
  [S, M, b] = StiffMassLoad(p, t, e, EqData, BdryData, time);

  u =                % to be completed!!!!!!!!!!!!!

  U=[U u];
  time = time + k/2;
  waitbar(time/FinalTime,h)

end


% visualization

figure(1)
set(gcf,'renderer','zbuffer')
h = plot(p,U(:,1));
set(gca,'XLim',[0 1],'YLim',[min(min(U)) max(max(U))])
for i=2:size(U,2)
  set(h,'YData',U(:,i))
  drawnow
end






% subroutines



function u0=InitialData(x)
  
  u0=zeros(size(x));



function [S, M, b] = StiffMassLoad(p, t, e, EqData, BdryData, time)

% Returns the global stiffness matrix S = (s_ij), mass matrix M = (m_ij)
% and load vector b = (b_j) with elements
% s_ji = int_Omega a Dphi_i Dphi_j, m_ji = int_Omega c phi_i phi_j and
% b_j = int_Omega f phi_j, respectively. 

% initializing

  n = size(p, 2);
  S = zeros(n, n);
  M = zeros(n, n);
  b = zeros(n, 1);
  S = sparse(S);
  M = sparse(M);

% collecting element contributions to S, M and b

  for el = 1:size(t, 2)
    [dS, dM, db] = ElementContributions(p, t, EqData, time, el);
    nn = t(1:2, el);
    S(nn, nn) = S(nn, nn) + dS;
    M(nn, nn) = M(nn, nn) + dM;
    b(nn)     = b(nn)     + db;
  end

% collecting boundary contributions to S, M and b

  for bel = 1:size(e, 2)
     [dS, db] = BoundaryContributions(p, e, BdryData, time, bel);
     nn = e(1, bel);
     S(nn, nn) = S(nn, nn) + dS;
     b(nn)     = b(nn)     + db;
  end




function [dS, dM, db] = ElementContributions(p, t, EqData, time, el)

% returns the element stiffness matrix dA and element load
% vector db for element number el defined by t(1:3,el) and
% p(t(1:2,el), and the equation data returned by EqData,
% using element midpoint quadrature.

  NodeCoords = p(t(1:2, el));
  x = sum(NodeCoords)/2;             % element midpoint!
  [a, c, f] = EqData(x, time, t(3, el));
  dx = diff(NodeCoords);             % element length
  phi = [.5 .5];
  Dphi = [-1/dx 1/dx];
  dS = [a*Dphi(1)*Dphi(1)*dx  a*Dphi(2)*Dphi(1)*dx;...
        a*Dphi(1)*Dphi(2)*dx  a*Dphi(2)*Dphi(2)*dx];
  dM = [c*phi(1)*phi(1)*dx  c*phi(2)*phi(1)*dx;...
        c*phi(1)*phi(2)*dx  c*phi(2)*phi(2)*dx];
  db = [f*phi(1)*dx; f*phi(2)*dx];



function [dA, db] = BoundaryContributions(p, e, BdryData, time, bel)

% returns the contributions dA and db from boundary element/node
% bel to the stiffness matrix and load vector.

  NodeCoord = p(e(1, bel));
  x = NodeCoord;
  [k, g] = BdryData(x, time, e(2, bel));
  phi = 1;
  dA = k*phi*phi;
  db = k*g*phi;



function [a, c, f] = EqData(x, t, tag)

% returns the values at x and time t of the data a=a(x,t), c=c(x,t)
% and f=f(x,t) in the differential equation cdu/dt-D(aDu)=f

  a = 1;
  c = 1;
  f = 0;



function [k, g] = BdryData(x, t, tag)

% returns the values at x of the data k=k(x,t) and g=g(x,t)

  if tag==2
     k = 10000; g = sin(5*t);
  elseif tag==1
     k = 10000; g = 0;
  end
  
