function [x,u]=newmarkwave(xspan,tspan,nstep,param,...
               c,u0,v0,g,f,varargin)
%NEWMARKWAVE risolve l'equazione delle onde
%       con il metodo di Newmark
%   [X,U]=NEWMARKWAVE(XSPAN,TSPAN,NSTEP,PARAM,C,...
%   U0,V0,G,F)
%   risolve l'equazione delle onde
%      D^2 U/DT^2 - C D^2U/DX^2 = F in
%   (XSPAN(1),XSPAN(2)) X (TSPAN(1),TSPAN(2)) con il
%   metodo di Newmark, con condizioni iniziali
%   U(X,0)=U0(X), DU/DX(X,0)=V0(X) e condizioni al
%   bordo di Dirichlet U(X,T)=G(X,T) in X=XSPAN(1)
%   ed in X=XSPAN(2). C e' una costante positiva,
%   F,G,U0 e V0 sono inline function.
%   NSTEP(1) e' il numero di intervalli di spazio.
%   NSTEP(2) e' il numero di intervalli in tempo.
%   PARAM(1)=ZETA e PARAM(2)=THETA.
%   [X,U]=NEWMARKWAVE(XSPAN,TSPAN,NSTEP,PARAM,C,...
%   U0,V0,G,F,P1,P2,...) passa i parametri addizio-
%   nali P1,P2,... alle funzioni U0,V0,G,F.
h  = (xspan(2)-xspan(1))/nstep(1);
dt = (tspan(2)-tspan(1))/nstep(2);
zeta = param(1);  theta = param(2);
N = nstep(1)+1;
e = ones(N,1); D = spdiags([e -2*e e],[-1,0,1],N,N);
I = speye(N);
lambda = dt/h;
A = I-c*lambda^2*zeta*D;
An = I+c*lambda^2*(0.5-zeta)*D;
A(1,:) = 0; A(1,1) = 1; A(N,:) = 0; A(N,N) = 1;
x = linspace(xspan(1),xspan(2),N);
x = x';
fn = feval(f,x,tspan(1),varargin{:});
un = feval(u0,x,varargin{:});
vn = feval(v0,x,varargin{:});
[L,U]=lu(A);
alpha = dt^2*zeta; beta = dt^2*(0.5-zeta);
theta1 = 1-theta;
for t = tspan(1)+dt:dt:tspan(2)
    fn1 = feval(f,x,t,varargin{:});
    rhs = An*un+dt*I*vn+alpha*fn1+beta*fn;
    temp = feval(g,[xspan(1),xspan(2)],t,varargin{:});
    rhs([1,N]) = temp;
    u = L\rhs;    u = U\u;
    v = vn + dt*((1-theta)*(c*D*un/h^2+fn)+...
        theta*(c*D*u/h^2+fn1));
    fn = fn1;    un = u;    vn = v;
end