function rSystem = odefun(t,x,d,z,p,St,Co,fnames)
%rSystem = odefun(t,x,d,z,p,St,Co,fnames)
%INPUTS: 
%t = current time
%x = current continous state
%d = current continous state derivative
%z = current system discrete state cell array
%p = parameter cell array
%St = state assignment to modules
%Co = connectivity matrix defnining shared states
%fnames = function handles of ode functions of system modules
%MODULES: 
%   - module 1: core network w voltage sensor and power sensor
%   - module 2: timer for tap changer
%   - module 3: timer for capacitor bank
%   - module 4: load
%   - module 5: capacitor bank
%   - module 6: wind turbine

%get the number of modules
nMo = size(St,1);
%get the number of states
nSt = length(x);
%initialize state cell array
xMod = cell(nMo,1);
%initialize state derivative cell array
dMod = cell(nMo,1);
%initialize full system residue
rSystem = zeros(nSt,1);

ni = 0; %running index for full system residues
for i = 1:nMo
    %collect current state and state derivative values for module i
    xMod{i} = x(St(i,1):St(i,2));
    dMod{i} = d(St(i,1):St(i,2));
    %compute residue for module i by calling that modules ode function
    rCur = fnames{i,1}(t,xMod{i},dMod{i},z{i},p{i});
    %compute the number of residues returned by module i
    nCur = length(rCur);
    %update full system residue 
    rSystem(ni+1:ni+nCur) = rCur;
    %update running index
    ni = ni + nCur;
end

%the residues of individual modules are no computed
%what remains is the residues corresponding to interconnection of modules
c = size(Co,1); %number of connections
for i = 1:c
    %residue  = difference between identical state appearing in different
    %modules
    rSystem(ni+i) = xMod{Co(i,1)}(Co(i,2)) - xMod{Co(i,3)}(Co(i,4));
end

