% 1d wave eqn:  u_tt = c u_xx
c=.02; n=100; h=1/n; dirichlet=1;
if dirichlet % Dirichlet BC
    x=[h:h:1-h]'; n1=n-1; 
    u0=exp(-100*(x-.5).^2);  % initial profile; initial velocity=0
    %u0=[x(1:75)/x(75); (1-x(76:end))/(1-x(75))];
    A=2*(1-c)*eye(n1)+c*diag(ones(n-2,1),1)+c*diag(ones(n-2,1),-1);
    u=(A/2)*u0; oldu=u0; plot(x,u0); axis([0 1 -1.2 1.2]); pause
    for i=1:2500
        newu=A*u-oldu; oldu=u; u=newu;
        if mod(i,10)==0 plot([0;x;1],[0;u;0]); axis([0 1 -1.2 1.2]); pause(.1); end
    end
else % neumann bc
    x=[0:h:1]'; n1=n+1;
    u0=exp(-100*(x-.5).^2);
    A=2*(1-c)*eye(n1)+c*diag(ones(n,1),1)+c*diag(ones(n,1),-1);
    A(1,2)=2*c; A(n1,n)=2*c;
    u=(A/2)*u0; oldu=u0; plot(x,u0); axis([0 1 -1.2 1.2]); pause
    for i=1:2500
        newu=A*u-oldu; oldu=u; u=newu;
        if mod(i,10)==0 plot(x,u); axis([0 1 -1.2 1.2]); pause(.1); end
    end
end