function [x,cal,b,exflg,time,Uout] = fsamp(fun,U,x0,ncal)
%
% This function samples the space defined in the matrix A. The output of
% this sampling is given in the function x.
%
% The needed input for this function are
%   fun @(x,inp,opt):   Function in which the non-linear system of
%                       equations is solved. In addition, the reassigment
%                       of the values of the Matrix A has to be performed
%                       in this function. The output of this function must
%                       at least be given by
%                       [x,
%
%   U:                  A matrix defining the different input values. The
%                       different values for each input is defined in the
%                       columns, whereas the rows define the different
%                       inputs.
%
%   x0:                 Solution of the system around the nominal point
%
%   ncal:               Number of calculated values in the function fun.
%                       The calculated variables have to be defined in fun
%                       as global variables in a matrix Cal (1xncal matrix)
%
% The output is defined as
%   x:                  Array of dimension length(A) x nout defining the
%                       output values for the input values defined in
%                       matrix A.
%
%   cal:                Output values which are calculcated in the
%                       simulation. This mostly include states. It requires
%                       the definition of the global variable Calg in fun.
%
%   b:                  Residuals in the Newton-Raphson approach for each
%                       screening point.
%
%   exflg:              Exitflag of fsolve for each screening point.
%
%   d:                  Cell structure containing the calculated results in
%                       DataTable. The size is defined by length(x0).
%
%   time:               Time needed at each iteration step for solving

global Calg

% Reset of Calg
Calg = zeros(1,ncal);

% Calculation of lengths
nsamp = size(U,1);

% Preallocation of memory to the values
xcell = cell(nsamp,1);
Uoutcell = xcell;
Uincell = xcell;
timecell = xcell;
acell = xcell;
bcell = xcell;

% Definition of the matrix U into a Cell structure
Uaux = mat2cell(U,ones(1,nsamp),size(U,2));
Uincell{1} = Uaux{1};
for k = 2:nsamp
    Uincell{k} = cell2mat(Uaux(k-1:k));
end

% Start of the timing
tic;
[xcell{1},acell{1},bcell{1},Uoutcell{1}] = fun(x0',Uincell{1});
xcell{1} = xcell{1}';
timecell{1} = toc;

% Solving of the system using a single for loop
parfor k = 2:nsamp
% for k = 2:nsamp
    
    % Start of the timing
    in = tic;
        
    % Solving of the non-linear system
    [xcell{k},acell{k},bcell{k},Uoutcell{k}] = fun(x0',cell2mat(Uincell(k)));
    xcell{k} = xcell{k}';

    % Stop of the timing
    timecell{k} = toc(in);
end
Uout = cell2mat(Uoutcell);
time = cell2mat(timecell);
x = cell2mat(xcell);
b = acell;cell2mat(acell);
exflg = cell2mat(bcell);
cal = 0;