function z=si(t)
% SI computes the function Si(t), which is the integral of sin(x)/x 
% from x=0 to x=t. This is called the sine integral.

y=ones(size(t));
nn=find(t~=0);
y(nn)=sin(t(nn))./t(nn);
z=cumtrapz(t,y);