Simulating LTI Systems in Matlab

This tutorial will assume you have created state-space matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ which characterize the system modeled by

$$ \mathbf{\dot{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t) $$ $$ \mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}(t) $$

and have stored these matrices in MATLAB variables A, B, C, and D, respectively. The time-domain response of the system can then be solved numerically in MATLAB using two functions, ss() and lsim(). The first few lines of the help documentation for each is given below

 >> help ss
	SYS = ss(A,B,C,D) creates an object SYS representing the continuous-time state-space model
		  dx/dt = Ax(t) + Bu(t) 
		   y(t) = Cx(t) + Du(t)
			 
 >>help lsim
	lsim(SYS,U,T) plots the time response of the dynamic system SYS to the input signal described by U and T. 
	The time vector T is expressed in the time units of SYS and consists of regularly spaced time samples. 
	The matrix U has as many columns as inputs in SYS and its i-th row specifies the input value at time T(i). For example, 
		t = 0:0.01:5;   u = sin(t);   lsim(sys,u,t)  
	simulates the response of a single-input model SYS to the input  u(t)=sin(t) during 5 time units.
			 

For a fully LTI system, where all state-space matrices are constant, the solution at time $t>0$ can be solved numerically by simply integrating $\mathbf{\dot{x}}(t)$ over time, beginning from some initial condition. For example, for the problem of homework #2, one way to setup the simulation is as follows

 
sys = ss(A,B,C,D); 				%combines matrices to generate a system 

nPeriods = 200;					%number of periods to simulate
pts_per_period = 10000;				%number of time-points per period
npts = nPeriods*pts_per_period;			%total number of time-points in simulation
t = linspace(0, nPeriods*Ts, npts);		%time vector for simulation
vp_Ts = [-Vg*ones(1,pts_per_period/2),...
		Vg*ones(1,pts_per_period/2)];	%creates one period of squarewave vp(t)
vp = repmat(vp_Ts, 1, nPeriods);		%repeats the single period of vp(t) nPeriods times
vs = vp/Vg*V;					%creates vs(t).  This just re-uses vp(t) and changes the amplitude 
vp = circshift(vp, round(tphi/Ts*pts_per_period)); 	%adds a phase-shift of tphi to vp

u = cat(1,vp,vs);				%concatenate vp and vs together for the u(t) vector needed in lsim
[Y,T,X] = lsim(sys,u,t,zeros(nstates,1));	%simulates the system with the input u(t), beginning at zero initial conditions