function [t,varargout]=rk4_multidim(f,dim,tinit,xinit,tfin,nstep)
%
% euler_multidim: 
% subroutine for a very naive implementation of the fourth order
% Runge-Kutta 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 of a row vector.
%             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 the time column vector.
%                                         % Horizontally concatenate the
x=cat(2,xinit,zeros(dim,nstep));          % the xinit row column vector  
%                                         % with the matrix where the time
%                                         % evolution will be allocated.
args1={};                                 % Initialize cell.
args2={};                                 % Initialize cell.
for i=1:dim,
    args1=cat(1,args1,xinit(i));          % Data allocation into cell.
    args2=cat(1,args2,xinit(i));          % Data allocation into cell.
end
%                                         % Start of the main loop.  
for i=1:nstep,
    t(i+1)=t(i)+h;                        % Update time.
    rk1=h*feval(f,t(i),args1{:});         % Pass args to f.  
    for j=1:dim,
        args2{j}=args1{j}+rk1(j)/2;       % Update cell.
    end
    rk2=h*feval(f,t(i)+h/2,args2{:});     % Pass args to f.
    for j=1:dim,
        args2{j}=args1{j}+rk2(j)/2;
    end
    rk3=h*f(t(i)+h/2,args2{:});           % Pass args to f.
    for j=1:dim,
        args2{j}=args1{j}+rk3(j);         % Update cell.
    end
    rk4=h*f(t(i)+h,args2{:});             % Pass args to f.   
    x(:,i+1)=x(:,i)+(1/3)*(rk1/2+rk2+rk3+rk4/2);    % Update solution.
   for j=1:dim,
        args1{j}=x(j,i+1);                % Update cell.
   end
end
%                                         % Main loop is over.
    varargout{1:dim,:}=x(1:dim,:);        % Store solution in varargout.
end
