function [uT]=solv_1d_wave(u0, ut0, T)
%
N=size(u0,2)-1;
u0_odd=[u0(1:N), -u0(N+1:-1:2)];
cf_u0=imag(fft(u0_odd));
    % Fourier coefficients in sine expansion of u0
ut0_odd=[ut0(1:N), -ut0(N+1:-1:2)];
cf_ut0=imag(fft(ut0_odd));
    % Fourier coefficients in sine expansion of ut0
cf_uT=cf_u0(2:N+1).*cos([1:N]*pi*T)+cf_ut0(2:N+1).*sin([1:N]*pi*T)./[1:N]*pi;
    % Fourier coefficients in sine expansion of u(x,T)
uT=real(ifft([0, sqrt(-1)*cf_uT, -sqrt(-1)*cf_uT(N-1:-1:1)]));
uT=uT(1:N+1);
%
%
