function [sol] = int_sys(func,var0,t,sub)

% This function integrates the DAE defined previously.
% The output of this function is defined in the structure sol in the
% following way:
%   sol.x:   	matrix in which the rows corresponds to the differential 
%               states and the columns corresponds to the time
%   sol.z:      matrix in which the rows corresponds to the algebraic 
%               states and the columns corresponds to the time
%   sol.J:      number corresponding to the overall cost function
%
% The required input of this function is given as
%   func:       Function giving the output of an integrator step
%   var0:       Structure containing:
%               x   initial values for the differential states
%               z   initial values for the algebraic states
%               p   values for the input
%                   if p is a vector, the same input is applied over the
%                   whole integration domain
%                   if p is a matrix, it needs to have t.N columns
%                   corresponding to the different time steps
%   t:          Structure of time parameters, containing 
%               int time interval inbetween the steps
%               N   total number of steps
%   sub:        Values for substitution

% Definition of the lenght of the systems
n_x = length(var0.x);
n_z = length(var0.z);

% Import of Casadi
import casadi.*

% Predefinition of the outlet vectors
sol.x = zeros(n_x,t.N+1);   sol.x(:,1) = full(var0.x);
sol.z = zeros(n_z,t.N+1);   sol.z(:,1) = full(var0.z);
sol.J = 0;


if size(var0.p,2) == 1
    sol.p = var0.p*ones(1,t.N+1);
    sol.p(4:end,:) = sol.p(4:end,:)+sub;
else
    sol.p = var0.p;
end

sub_size = size(sub,2);

% Calculation of the new values
for j = 2:t.N+1
    
    if sub_size > 1
        if sub(1,j) ~= sub(1,j-1)
            fprintf('Change of flowrate difference from %3i to %3i after %3i min. \n', sub(1,j-1),sub(1,j),(j-1)*t.int/60);
        elseif sub(2,j) ~= sub(2,j-1)
            fprintf('Change of pressure difference from %3i to %3i after %3i min. \n', sub(2,j-1),sub(2,j),(j-1)*t.int/60);
        elseif sub(3,j) ~= sub(3,j-1)
            fprintf('Change of temperature difference from %3i to %3i after %3i min. \n', sub(3,j-1),sub(3,j),(j-1)*t.int/60);
        elseif sub(4,j) ~= sub(4,j-1)
            fprintf('Change of composition difference from %2i wt.%% to %2i wt.%% after %3i min. \n', sub(4,j-1)*100,sub(4,j)*100,(j-1)*t.int/60);
        end
    end
    func_eval = func('x0',sol.x(:,j-1),'z0',sol.z(:,j-1),'p',sol.p(:,j-1));
    sol.x(:,j) = full(func_eval.xf);
    sol.z(:,j) = full(func_eval.zf);
    sol.J = sol.J + full(func_eval.qf);
end

fprintf('The dynamic system is modelled.\n')
