alpha=0.5;
N=1000;  % must be even!
x=linspace(0,1,N);
ysteady=(1+1/2/alpha/alpha)*x-x.*x/2/alpha/alpha;
fftysteady=fft(ysteady);
k=[0:N/2 -(N/2-1):-1]; %N frequencies
ytransient=-real(ifft(exp(-alpha*alpha*k.*k*t).*fftysteady));
yt=ysteady+ytransient;
plot(x,yt)