function [y,x] = rc_sim(mttx_0,mttU,mttT,MTTpar);
## Initialise
[MTTu] = zero_input(1); # Zero the input
[MTTx] = mttx_0; # Read in initial state
N = length(mttT);
Nx = length(mttx_0);
y = zeros(1,N);
x = zeros(Nx,N);
mttdt = mttT(2)-mttT(1);
for i = 1:N
MTTt = mttT(i);
MTTu = mttU(i);
[MTTy] = rc_cseo(MTTx,MTTu,MTTt,MTTpar); # Output
y(1,i) = MTTy;
x(:,i) = MTTx;
if 1>0 # Dont if no states
[MTTdx] = rc_cse(MTTx,MTTu,MTTt,MTTpar); # State derivative
[mttAA] = rc_smxa(MTTx,MTTu,mttdt,MTTpar); # (I-Adt)
[mttAA] = reshape(mttAA,Nx,Nx);
[mttAAx] = rc_smxax(MTTx,MTTu,mttdt,MTTpar); # (I-Adt)x
[MTTopen] = rc_switchopen(MTTx); # Open switches
[MTTx] = mtt_implicit(MTTx,MTTdx,mttAA,mttAAx,mttdt,1,MTTopen); # Implicit update
else
endif; # 1>0
endfor; # Integration loop
endfunction