function U=NewmarkFEM(f,g,u0,udot0,T,c,k,Nt)

% U=NewmarkFEM(f,g,u0,udot0,T,c,k,Nt)
%
% Input:
%     f    : source term (vectorized function of x,y,t)
%     g    : Neumann condition (vector valued, vectorized)
%     u0   : initial condition (function of x,y)
%     udot0: initial condition for u_t (function of x,y)
%     T    : triangulation
%     c    : speed of waves in the medium
%     k    : time-step
%     Nt   : number of time-steps
% Output:
%     U    : U(i-1,n-1) : solution at x_i,t_n

[S,M,Mlump,T,Dir,Free]=FEMprepare2wave(T);
S=c^2*S(Free,Free);
M=M(Free,Free);
U=zeros(size(T.coordinates,1),Nt+1);
x=T.coordinates(:,1);
y=T.coordinates(:,2);
U(:,1)=u0(x,y);
V0=udot0(x,y);
U(Free,2)=U(Free,1)+k*V0(Free)...
          +k^2/2*(M\(loadplustraction(T,f,g,0,Free)-S*U(Free,1)));

for n=2:Nt     % n=current time
    U(Free,n+1)=2*U(Free,n)-U(Free,n-1)...
                +k^2*(M\(loadplustraction(T,f,g,(n-1)*k,Free)-S*U(Free,n)));
end

retur