¿Cómo resolver un circuito de refuerzo (o cualquier circuito RLC conmutado) en el dominio del tiempo?

Necesito modelar un circuito boost con la siguiente topología:

^{simular este circuito : esquema creado con CircuitLab}

Estoy interesado en el voltaje de la resistencia, así que encontré la ecuación diferencial de cada estado y los resolví por el método implícito de Euler, alternando la ecuación a resolver a medida que cambiaba el estado.

Comparé los resultados con un simulador de circuito. Cada ecuación diferencial es correcta por sí sola, pero cuando las pongo juntas, una tras otra, formando la dinámica conmutada de refuerzo, el resultado es incorrecto.

Estos son mis pasos:

Ecuación diferencial de 'estado ON' (solo la descarga del capacitor)

$$\frac{dV_{out}}{dt} = -\frac{V_{out}}{RC}$$

Ecuación diferencial de 'estado APAGADO'

$$\frac{d^2V_{out}}{dt^2} + \frac{1}{RC}\frac{dV_{out}}{dt} + \frac{1}{LC}V_{out} = \frac{V_{in}}{LC}$$

Este es mi código (MATLAB):

R = 10; 
L = 1e-5; 
C = 1e-4;
Vin = 10; 

h = 1e-6;   % time step
sizeVector = 10000; % number of steps
t = 0:h:sizeVector*h-h;  % time vector

Vout = zeros(1,sizeVector); % output voltage vector
dVout = zeros(1,sizeVector); % first derivative vector

freqSw = 5e3; % switching frequency
D = zeros(1,sizeVector); % State of switch for each simulation step (duty cycle = 0.5)
D = sign(cos(2*pi*freqSw.*t + pi/2 + 0.001));
D(D<0) = 0;

for k=2:sizeVector
    if D(k) == 0
        % Off state
        Vout(k) = (Vout(k-1) + h*((dVout(k-1) + (h/(L*C))*Vin)/(1 + h/(R*C))))/(1 + (h^2/(L*C))/(1 + h/(R*C)));
        dVout(k) = (dVout(k-1) + h*(Vin/(L*C) - Vout(k)/(L*C)))/(1 + h/(R*C));
    end
    if D(k) == 1
        % On state
        Vout(k) = Vout(k-1)/(1+h/(R*C));
        dVout(k) = (Vout(k)-Vout(k-1))/h; % just for convenience
    end

end


figure(1)
plot(t,Vout)
xlabel('Time (s)')
ylabel('Voltage (V)')
hold on

Estas son las comparaciones entre los resultados de PSIM (rojo) y mi código en MATLAB (azul). (Para la descarga del capacitor consideré en ambos casos una tensión inicial de 10V)

Estado abierto:

Estado cerrado:

Boost convertidor como en el código anterior:

Nunca puedo obtener un voltaje de estado estable superior a mi voltaje de entrada (10 V, en este caso), pero la forma de onda parece ser consistente con cada una de las formas de onda de los estados por separado.

¿Debe haber un error en mi código o me estoy perdiendo algún punto importante? ¿Es posible este tipo de resolución?

Editado después de las respuestas:

Siguiendo las sugerencias, cambié la ecuación diferencial de segundo orden en un sistema de dos ecuaciones de primer orden, para calcular la corriente del inductor, lo cual es necesario porque cuando comienza el 'estado apagado', la corriente del inductor debe ser una condición inicial. En el modelo anterior, no estaba considerando el cambio actual en el 'estado activado'.

El sistema 'fuera del estado':

$$\begin{equation} \begin{bmatrix} \frac{dV_{out}}{dt}\\ \frac{dI_{ind}}{dt} \end{bmatrix} = \begin{bmatrix} -\frac{1}{RC} & \frac{1}{C}\\ -\frac{1}{L} & 0 \end{bmatrix} \times \begin{bmatrix} V_{out}\\ I_{ind} \end{bmatrix} + \begin{bmatrix} 0\\ \frac{Vin}{L} \end{bmatrix} \end{equation}$$

El sistema 'en estado': (en realidad no es un sistema, pero lo mantengo matricial por simplicidad)

$$\begin{equation} \begin{bmatrix} \frac{dV_{out}}{dt}\\ \frac{dI_{ind}}{dt} \end{bmatrix} = \begin{bmatrix} -\frac{1}{RC} & 0\\ 0 & 0 \end{bmatrix} \times \begin{bmatrix} V_{out}\\ I_{ind} \end{bmatrix} + \begin{bmatrix} 0\\ \frac{Vin}{L} \end{bmatrix} \end{equation}$$

Todavía hay un problema. De esta manera, no puedo evitar que la corriente del inductor sea negativa en el 'estado apagado', por lo que inserté una declaración if en cada iteración del 'estado apagado'. Si$I_{ind} < 0$, resolver el sistema nuevamente imponiendo$I_{ind} = 0$. Esto se puede hacer simplemente poniendo a cero los elementos relacionados con$I_{ind}$en la matriz.

A continuación se muestra la comparación de resultados (PSIM, primero, MATLAB segundo) y el código.

R = 10; 
L = 1e-3; 
C = 1e-4;
Vin = 10; % DC input voltage

h = 1e-7;                % step
sizeVector = 120000;     % number of steps
t = 0:h:sizeVector*h-h;  % time vector

X = zeros(2,sizeVector); % state vector 

A = [-1/(R*C) 1/C ; -1/L 0]; % off state matrix
b = [0; Vin/L];              % off state vector

Aaux = [-1/(R*C) 0 ; 0 0];   % off state auxiliar matrix
baux = [0 ; 0];              % off state auxiliar vector

A2 = [-1/(R*C) 0; 0 0];      % on state matrix
b2 = [0; Vin/L];             % on state vector

freqSw = 5e3;               % switching frequency
D = zeros(1,sizeVector); % State of switch for each simulation step (duty cycle = 0.5)
D = sign(cos(2*pi*freqSw.*t + pi/2 + 0.001));
D(D<0) = 0;

for k=2:sizeVector

    if D(k) == 0 % off state 
        X(:,k) = (inv(eye(length(A)) - h*A))*(X(:,k-1) + h*b);
        if X(2,k) < 0 % avoid Iind < 0
            X(:,k) = (inv(eye(length(Aaux)) - h*Aaux))*(X(:,k-1) + h*baux);
        end

    end
    if D(k) == 1 % on state
        X(:,k) = (inv(eye(length(A2)) - h*A2))*(X(:,k-1) + h*b2);
    end
end


figure(1)
plot(t,X(1,:))
hold on
plot(t,X(2,:))
xlabel('Time (s)')
ylabel('Voltage (V)\Current(A)')
legend('V_{out}','I_{in}')
ylim([-6 36])
grid
hold on