N=100;
L=1;
u=1; %wave speed
x=linspace(0,L,N+1);
y0=exp(-50*(x/L-0.5).^2); %initial wave shape
y0(N+1)=[]; % redundant value at endpoint
%y0=zeros(1,N);
%v0=exp(-100*(x/L-0.5).^2); %initial velocity pulse shape
%v0(N+1)=[]; % redundant value at endpoint
v0=zeros(1,N);

y0ext=[0 y0(2:N) 0 -y0(N:-1:2)]; %extend to force a sine series
v0ext=[0 v0(2:N) 0 -v0(N:-1:2)]; %extend to force a sine series
k=pi/L*[0:N (-N+1):-1]; %frequencies
omega=u*k; % note omega(1)=0:  makes slight problem below;
A=fft(y0ext); %Fourier coefficients of sine series of y0
omega(1)=1; % temporarily!  To avoid a divide by zero error in the line below.
B=fft(v0ext)./omega; %B(1) is actually zero, so changing omega(1) is irrelevant
omega(1)=0; % undo change
yext=real(ifft(A.*(cos(omega*t))+B.*sin(omega*t))); %ifft to get y, and enforce "real"
plot(x,yext(1:N+1)) % just plot the part corresponding to original interval