Formulate a set of state equations with the form:
dX/dt = [A]X(t)+B*vin(t)
where X is a vector with 5 elements (Vc1,Vc2,Vc2,Il1,Il2), [A] is a
5x5 real matrix, and B is a vector with 5 elements. This is quite easy.
Solve it by numerical integration with the trapezoidal rule:
X(t0+dt)=X(t0)+(dt/2)*([A]*X(t0)+B*vin(t0)+[A]*X(t0+dt)+B*vin(t0+dt))
where t0 is the present time, dt is a small time interval ahead, X(t0)
is the present state, and X(t0+dt) is the next state.
The required calculation is:
X(t0+dt)=[[I]-(dt/2)*[A]]^(-1)*(dt/2)*([A]*X(t0)+B*(vin(t0)+vin(t0+dt)))
The matrix inversion has to be done just once.
This results in good precision if you put about 50 "dt" for each cycle.
If I find time I will make a starting code with these calculations