Sum function handles efficiently

I need to solve an ODE where the motion is determined by charges in position (where and ). For the sake of simplicity let's just assume that and that the dynamic given by the particle in is given by (i.e. I need to solve). I just need to solve an ODE where the motion is determined by the effects of the charges combined, and in my case the dynamic is simply given by the sum of the functions . I tried definining function handles stored in a 2 dimensional array of size , where in position I store the function . Then, I would need to define and solve an ODE where the motion is determined by F. The way I did this is by recursion, i.e.:

%I have already defined f as a 2D array, where f{i,j}=f_{i,j} described
%in the text
F= @(t,x) 0;
for i= 1:1:N
    for j=1:1:N
        F= @(t,x) F(t,x) + f{i,j}(t,x)
    end
end

After this, I use the solve function:

fun = ode(ODEFcn=@(t,x) F(t,x),InitialTime=0,InitialValue=[0,0]);    % Set up the problem by creating an ode object
sol = solve(fun,0,100);                                               % Solve it over the interval [0,10]

The problem is: the performance is very bad. I already see this when defining F. I think there might be some issues with the recursion, nad maybe there's a better way for defining F, in such a way that the performances get better.

[Torsten]

Are the f{i,j} complicated ?

If not, don't define them as separate N^2 function handles, but just compute their values in a normal function F and sum them.

[Stephen23]

"The problem is: the performance is very bad."

That is not surprising: you are creating lots of separate workspaces and calling lots of functions. Ouch.

"... maybe there's a better way for defining F, in such a way that the performances get better."

Why not simply write a normal function that calls f in a loop (or two) and sums the outputs?:

https://en.wikipedia.org/wiki/KISS_principle

[Tommaso Cortopassi]

@Stephen23 @Torsten yeah, I'll do that. I tried before using a function with a for cycle but for some reasons I got an error in "solve", for not passing enough arguments. Now I tried again and it seems to work. I don't know, for some reason I thought that having a long complicated function handle would have been computationally more efficient than a function with a for cycle, even though I am probably mistaken

[Alex]

Not sure if this makes it any better, but you can avoid the recursion and define your combined ODE in two lines:

[i,j] = meshgrid(1:N, 1:N);
F = @(t,x)sum(arrayfun(@(i,j)f{i,j}(t,x),i(:),j(:)));