function [u2]=newmark(N,M,tau)
%
%  Schema de Newmark pour l'equation des ondes
%
%  parametres
%
%  N        : nombre de points interieurs dans l'intervalle [0,1]
%  h        : pas d'espace
%  M        : nombre de pas de temps
%  tau      : pas de temps
%  t        : temps courant
%  u0     : N-vecteur, u0(i) est une approximation de u(x_i,t_n-1)
%  u1     : N-vecteur, u1(i) est une approximation de u(x_i,t_n)
%  u2     : N-vecteur, u2(i) est une approximation de u(x_i,t_n+1)
%
h=1./(N+1);
lambda=tau^2/h^2;
%
% condition initiale u0 et u1
%
for i=1:N
  u0(i)=???;
end
u1(1)=???;
for i=2:N-1
  u1(i)=???;
end
u1(N)=???;
%
% schema de Newmark
%
t=tau;
for n=2:M
  t=t+tau;
  u2(1)=???;
  for i=2:N-1
    u2(i)=???;
  end
  u2(N)=???;
%
% reactualiser la solution et imprimer la norme euclidienne de u2
%
  norm2=0;
  for i=1:N
    u0(i)=u1(i);
    u1(i)=u2(i);
    norm2=norm2+u2(i)*u2(i);
  end
  fprintf(' pas de temps %d temps %f norm2 %e \n',n,t,norm2)
end
%
% imprimer l'erreur maximum au temps final
%
err = 0;
for i=1:N
    erri = abs(u2(i)-uex(i*h,t));
    if (erri>err) 
        err = erri;
    end
end
fprintf(' erreur maximum au temps final  %e \n',err)
%
% definition de la condition initiale w(x)
%
function init = w(x)
init = ???;    
%
% solution exacte (utilise pour le calcul de l'erreur)
%
function u_exact = uex(x,t)
u_exact =  0.5*(omega(x-t)+omega(x+t));
%
% definition de la condition initiale prolongee 
%          par imparite et 2-periodicite
%
function init = omega(x)
y = mod(x,2);
if (y<=1) 
  init = w(y);
else
  init = -w(y-1);
end
