t=0:.05:5; tau=1;
[nP1,dP1]=paderm(tau,1,3); GP1=tf(nP1,dP1)

[nP2,dP2]=paderm(tau,1,5); GP2=tf(nP2,dP2)

[nP3,dP3]=pade(tau,3); GP3=tf(nP3,dP3)

G=tf(1,[1,1]);   ytmp=step(G,t);
ii=find(t>tau); 
y0=[zeros(ii(1)-1,1); ytmp(1:length(ii))];
y1=step(GP1*G,t); y2=step(GP2*G,t); y3=step(GP3*G,t);
plot(t,y0,t,y1,'--',t,y2,':',t,y3,'-.')

set(G,'Td',tau); %set the time delay to the original system
y0=step(G,t);