function [t,varargout]=euler_multidim(f,dim,tinit,xinit,tfin,nstep)

% euler_multidim: 
% subroutine for a very naive implementation of the
% Euler integrtaion scheme for arbitrary dimension "dim"
%
% INPUT: 
% f:           name of a row vector field assigned in the form
%                        f(t,y1,...,ydim) where
%              "t" is time (the explicit form of f may not depend upon t)
%              "y1,...,ydim" are "dim" spatial coordinates.
%              f may be assigned as a "function handle" or as an
%              anonymous function.
% dim:         number of dimensions.
% tinit:       initial time.
% xinit:       initial position.
% tfin:        final time.
% nstep:       number of integration steps.
%
% OUTPUT
% t:           time.
% varargout:   MATLAB pre-defined statement. It is used only inside 
%              a function M-file to contain the optional output arguments 
%              returned by the function.
%
%
h=(tfin-tinit)/nstep;                     % Define the time mesh.
t=[tinit zeros(1,nstep)];                 % Define time column vector.
x=cat(2,xinit,zeros(dim,nstep));          % Concatenate horizontally
%                                         % the "xinit" row vector
%                                         % with dim*nstep rectangular  
%                                         % where the time evolution data
%                                         % will be stored.
%
args={};                                  % Initialize a cell.
%                                         % Store by vertical concatenation
for i=1:dim,                              % into the cell the values
    args=cat(1,args,xinit(i));            % of the initial data. Only data
end                                       % separated by a commas (cells)
%                                         % can be passed to f. 
%                                         
for i=1:nstep,                            % Main cycle
    t(i+1)=t(i)+h;                                  % Time update.
    x(:,i+1)=x(:,i)+h*feval(f,t(i),args{:});        % Solution update.
   for j=1:dim,                           % Update the cell: only 
        args{j}=x(j,i+1);                 % cell-ordered data can be  
   end                                    % passed to f.
end                                       % End of the main cycle.
%
%                                         % Store solution data into
    varargout{1:dim,:}=x(1:dim,:);        % varargout statement. 
end
