function y=wave1(h,k)

     lam=(k/h)^2;
     N=1/h+1;
     A=zeros(N);

% Dirichlet BCs
%     A(1,1)=1;
%     A(N,N)=1;


% Neumann BCs
     A(1,1)=(2-2*lam);
     A(1,2)=2*lam;
     A(N,N)=(2-2*lam);
     A(N,N-1)=2*lam;
% End of BCs 

     for i=2:(N-1)
        A(i,i)=2-2*lam;
        A(i,i-1)=lam;
        A(i,i+1)=lam;
     end

     y=zeros(N,100);


% Smooth IC
%for i=1:N
%     y(i,1)=exp(-((i-N/2)*h)^2/.001);
%end

% Discontinuous IC
for i=(N/2-100):(N/2+100)
     y(i,1)=1;
end

% Main Loop
for i=2:N-1
y(i,2)=(lam*y(i+1,1)+(2-2*lam)*y(i,1)+lam*y(i-1,1))/2+k*0;
end
y(1,2)=0;
y(N,2)=0;

     for i=2:99
y(:,i+1)=A*y(:,i)-y(:,i-1);;
end
end
