function[y,t]=simulate(a,b,c,u,t,x0,base,M,C,K,f)

% This function simulates the state space system given by
% dx/dt = a x + b u
%     y = c x
% at each time in t.  x0 refers to the initial conditions.
% base should be set to 1 for base excitation, 0 otherwise
% M,C, and K need only be specified for base exitiation


for i=1:length(t)
  xh(i,:)=(expm(a*t(i))*x0')';
end

[m,n]=size(a);
[m,nb]=size(b);
s=expm([[a b]*(t(2)-t(1)); zeros(nb,n+nb)]);
ad=s(1:n,1:n);
bd=s(1:n,n+1:n+nb);


xp=zeros(1,length(a));
for i=1:length(t)-1
  xp(i+1,:)=(ad*xp(i,:).' + bd*u(i)).';
end

x=xh+xp;

if base
   [N,N]=size(M);
  

  for i=1:length(t)
     xa(i,1:N)=(inv(M)*(f'*u(i)-C*x(i,N+1:2*N)'-K*x(i,1:N)')+u(i))';
     %x=xa;  
  end
  x=[xa xa];
end
y=(c*x')';
