matrix numerical-methods ode scilab state-space

How to solve a set of ODEs in Scilab?

I have following set of ODEs which I would like to numerically solve in Scilab

I have successfully written the function for evaluating the right hand side of the first equation in other words I am able to solve the first equation

function ut = u(t)
    ut = [Vm*cos(2*%pi*fs*t); Vm*sin(2*%pi*fs*t)];
endfunction

function dxdt = SystemModel(t, x, u)
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    dxdt = A*x + B*u(t);
endfunction

My problem is that I don't know how to write similar function for evaluation of the right hand side of the second equation because it depends on solution of the first equation. Can anybody give me an advice how to do that?

Possible solution:

x0 = zeros(4, 1);
xtilde0 = zeros(4, 1);
X0 = [x0; xtilde0];
t0 = 0;
dt = 0.001;
t = 0:dt:1;

function ut = u(t)
    ut = [Vm*cos(2*%pi*fs*t); Vm*sin(2*%pi*fs*t)];
endfunction

function dXdt = RightHandSide(t, X, u)

    x      = X(1:4);
    xtilde = X(5:8);

    // dx/dt = A*x + B*u 
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    // dxtilde/dt = (An - L*Cn)*xtilde + (dA - L*dC)*x + dB*u
    An = [(-RSn*(LLn + LMn)^2 - RRn*LMn^2)/(LLn*LMn*(LLn + LMn)),                                                     0.0, RRn/(LLn*(LLn + LMn)),             pp*wm/LLn;
                                                             0.0,  (-RSn*(LLn + LMn)^2 - RRn*LMn^2)/(LLn*LMn*(LLn + LMn)),            -pp*wm/LLn, RRn/(LLn*(LLn + LMn));
                                           (LMn*RRn)/(LLn + LMn),                                                     0.0,      -RRn/(LLn + LMn),                -pp*wm;
                                                             0.0,                                   (LMn*RRn)/(LLn + LMn),                 pp*wm,     -RRn/(LLn + LMn)];
                                                             
    K = 1.5;
    l1 = (K - 1.0)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LLn*LMn*(LLn + LMn)) + RRn/(LLn + LMn));
    l2 = (K - 1.0)*pp*wm;
    l3 = (K^2 - 1.0)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LMn*(LLn + LMn)) - (LMn*RRn)/(LLn + LMn)) - (K - 1)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LMn*(LLn + LMn)) + (LLn*RRn)/(LLn + LMn));
    l4 = -(K - 1.0)*LLn*wm*pp;
    L = [l1, l2; 
        -l2, l1; 
         l3, l4;
        -l4, l3];
        
    Bn = [(LLn + LMn)/(LLn*LMn), 0.0; 0.0, (LLn + LMn)/(LLn*LMn); 0.0, 0.0; 0.0, 0.0];
                                                             
    Cn = [1.0, 0.0, 0.0, 0.0; 0.0, 1.0, 0.0, 0.0];
    
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    C = [1.0, 0.0, 0.0, 0.0; 0.0, 1.0, 0.0, 0.0];
    
    dA = An - A;
    dB = Bn - B;
    dC = Cn - C;
    
    dxdt      = A*x + B*u(t);
    dxtildedt = (An - L*Cn)*xtilde + (dA - L*dC)*x + dB*u(t);
    dXdt = [dxdt; dxtildedt];
endfunction

X = ode(X0, t0, t, list(RightHandSide, u));

Solution

Let y = x_tilde. Let assume that it is a 3x1 vector (we can't guess its size with your current presentation).

Build the column X = [x1 x2 x3 x4 y1 y2 y3].' (big X)
Express the column (dX/dt) according to X coordinates and t
Convert the system built in 2) into a Scilab function X_dot = Xder(t, X)
Build the initial state vector Xinit = [x1(t_init); x2(t_init); .. y3(t_init)]
Define the vector t of times to which you want the values of X. They all likely have to be ≥ t_init, and must be strictly increasing.
call X = ode(Xinit, t_init, t, Xder)

X(:,i) should then be the values of X components for each t(i) date.

You can "back-split" big X into x = X(1:4,:) and x_tilde = X(5:$,:).