% Refer to "LiveLink TM for MATLAB ®. User’s Guide", section "Function Input/Output Considerations" :
% 1. When you write your own functions, remember that the input arguments are vectors.
% 2. Element-by-element multiplication, division, and power must be used—that is, the
% operators .* , ./ , and .^ .
% 3. Attention to "ifs" !!

% Also, If the function M-file is modified using a text editor, in COMSOL click "Clear Functions" to ensure that	
% the functions’ modifications are updated in the COMSOL Multiphysics model.

% Euler-Maruyama method for linear SDE
function Xt_calc =  EM_linearSystem(name, f, g, st, X0, T, N, dtMult, dimSystemEq, legend_str)

    % SDE: dX = f(Xt)  dt + g(Xt) dW
    % where both f and g are linear 

    if (nargin < 10)
        
        fprintf('ERROR: Not enough arguments');
        return;
    end
    
    % constants
    dt = T/N;
    Dt = dtMult*dt; L = N/dtMult;

    % initializations
    Xt = zeros(1,L);
    % do not randn('state',100): use rng() instead.
    % https://es.mathworks.com/help/matlab/math/updating-your-random-number-generator-syntax.html
    rng('default')
    dW = sqrt(dt)*randn(dimSystemEq,N);	% Brownian increments for each path
    % dW = normrnd(0,dt);  % TODO decide
    % dW = random(makedist('Normal', 0, sqrt(dt)));
    W = cumsum(dW);		% discretized Brownian path/s

    % main
    clear Xprev; clear Winc; clear Xt;
    
    Xprev = X0;     % dimSystemEq vector
    % Xprev = X0';  % dimSystemEq column vector
    for j = 1:L
        
        for WienerPath = 1:dimSystemEq  
            % Calculates each of the path/s, one per each dimension dimSystemEq
            Winc(WienerPath) = sum(dW(WienerPath, dtMult*(j-1)+1:dtMult*j));
        end
        
        Xprev = Xprev + ...
                f(st.c1X(Xprev), st.c2) * Dt + ...
                g(st.sigma1X(Xprev), st.sigma2) * Winc';
            
        Xt(:, j) = Xprev;  % Xt(dimSystemEq, j)
    end
    
    % plots
    h = figure ('Name', name);
    ylim([0, inf]);
    for eqN = 1:dimSystemEq
        
        plot([0:Dt:T],[X0(eqN),Xt(eqN,:)]);
        xlabel('t','FontSize',12)
        ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
        hold on
    end
    
    legend( legend_str ); 
    hold off

    saveas(h, name);
    
    Xt_calc = Xt;
end