function aveImplicit
L = 3;   % length of the space interval
n = 350; % number of interior grid points
r = 0.8; % ratio to compute time step
%r = 1.8;% ratio to compute time step
T = 5;   % final time

function r = iv(x);  % initial displacement
  r = min([2*x,2-2*x]')';  r = max([r,0*x]')';
end%function

dx = L/(n+1);
dt = r*dx;
x  = linspace(0,L,n+2)';

An = spdiags([-ones(n,1), 2*ones(n,1), -ones(n,1)],[-1,0,1],n,n)/dx^2;
M0 = speye(n)+dt^2/4*An;
M1 = 2*speye(n)-dt^2/2*An;

y0 = iv(x);  y0 = y0(2:n+1);
y1 = y0;   % use zero initial speed
y2 = y0;
for t = 0:dt:T+dt;
  plot(x,[0;y0;0]); axis([0,L,-1,1]); drawnow()
  y2 = M0\(M1*y1-M0*y0);
  y0 = y1; y1 = y2;
end%for

end%function
