I would use one of the ode functions, e.g. ode45. You can then define an arbitrarily complex, nonlinear, time-varying, etc. system.
In your case, something like that (parameters are rubbish, but you have those)
function [Tout, Xout] = run_ode
% Some Parameters
Rm = 1000;
Ra = 1000;
% Time-varying parts:
C = @(t)(sqrt(t+1)^3);
dC = @(t)(3*sqrt(t+1)/2);
r = @(x)x*(-x<0); % Ramp function
p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];
% Set up simulation
dx = @(t,x)ode_sys(t,x,C(t),dC(t),p(x));
Ts = 1e-12;
t_end = 1e-10;
t_sim = 0:Ts:t_end;
x0 = randn(5,1);
[ Tout, Xout ] = ode45(dx,t_sim,x0);
end
% Define system equations
function dy = ode_sys(t,x,C,dC,p)
% Some Parameters
Rs = 100;
Cr = 1e-6;
Cs = 2*Cr;
Ca = 0.1*Cr;
Ls = 6e-9;
Rc = 1e6;
b1 = 1/(Rs*Cr);
b2 = 1/(Rs*Cs);
% Defined time-varying system matrices elementwise
Ac = [ -dC/C 0 0 0 0 ;
0 -b1 b1 0 0 ;
0 b2 -b2 0 1/Cs;
0 0 0 0 -1/Ca;
0 0 -1/Ls 1/Ls -Rc/Ls ];
Pc = [ 1/C -1/C; -1/Cr 0; 0 0 ; 0 1/Ca ; 0 0 ];
dy = Ac*x + Pc*p;
end
Cheers
