%
%  This routine determines the equivalent discrete-time transfer function
%  assuming a zero order hold
%
%  The routine then plots th response of a continuous-time plan and
%  it's discrete-tiem equivalent
%  
  clear

%
%  determine the final time
%
  Tf = 4.0;
%
%  determine the continuous-time time sequency and step input. For the
%  continuous-time system the time steps need to be closer together and do
%  not need to correlate with the discrete-time steps
%
  t = linspace(0,Tf,1000);
  u = ones(1,length(t));
%
%  determine the continuous-time plant
%
  Gp = tf(1,[1 3 2]);
%
%  simulate the continuos-time plant
%
  y = lsim(Gp,u,t);
%
%  set the sampling interval
%
  T = 0.25;
%
%  determine the discrete-time time sequence and step input
%
  tz = [0:T:Tf];
  uz = ones(1,length(tz));
 
  
%
%  convert to a discrete-time transfer function
%
  Gi = Gp * tf(1,[1 0]);
  Gp_z1 = c2d(Gi,T,'imp');  % use the impulse-invariant transformation
%
%  multiply by the delay transfer function and simplify
%
  Gpz = minreal(Gp_z1*tf([1 -1],[1 0],T));
%
%  determine the zeros of the transfer function
%
  zero(Gpz)
%
%  simulate the discrete-time system
%
  yz = lsim(Gpz,uz,tz);
%
%  make the staircase functions for plotting the discrete-time output
%
  [ts,ys] = stairs(tz,yz);
%
%  Now plot
%
  plot(t,y,'--',ts,ys,'-'); grid;
  legend('Continuous','Discrete');
  xlabel('Time (sec)');
  ylabel('y(t)');
