function w=wavefd(xl,xr,yb,yt,M,N)
% Ecuación de onda por diferencias avanzadas
%   Entrada: [xl,xr], tiempo[yb,yt], pasos M, tiempos N
 f=@(x) sin(pi*x); g=@(x) 0*x; l=@(t) 0*t; r=@(t) 0*t;
 c=2; close all
 h=(xr-xl)/M; k=(yt-yb)/N; m=M-1; n=N;
 sigma=c*k/h;
 a=diag(2-2*sigma^2*ones(m,1))+diag(sigma^2*ones(m-1,1),1);
 a=a+diag(sigma^2*ones(m-1,1),-1);                       % Matriz A
 lside=l(yb+(0:n)*k); rside=r(yb+(0:n)*k);
 w(:,1)=0.5*a*f(xl+(1:m)*h)'+k*g(xl+(1:m)*h)';           % Cond. iniciales
 w(:,2)=a*w(:,1)-w(:,1)+sigma^2*[lside(1); zeros(m-2,1); rside(1)];
 for j=2:n
    w(:,j+1)=a*w(:,j)-w(:,j-1)+sigma^2*[lside(j); zeros(m-2,1); rside(j)];
 end
 w=[lside; w; rside];                                    % + cond. contorno
 x=(0:m+1)*h; t=(0:n)*k;
 mesh(x,t,w');                                           % Plot 3-D   
 view(60,30); axis([xl xr yb yt -1 1])
end