Local function that makes anonymous functions

I recently discovered the joy of simulations and the wonder that is ode45. But it got me wondering how to make more dynamic systems as, least with my current methodology I have to write out the ODEs by hand, which for say a 100 particle electron cloud is impracticle. So, I gave myself the task of figuring out just this as a proof of concept, before going on to try some more involved simulations.

The best idea I could come up with is to make a function which spits out a collection of anonymous functions and starting conditions. I think I have it mostly made when I ran into an unexpected problem. It gave me this error.

Unrecognized function or variable 'x'.

I will give you a simplified ecample of what I am trying to do that made the error below as well as the full code in the comments for those interested. I am mostly looking for how to resolve the error, although if there is just a blatantly better way to go about this whole idea that is also welcome. (Honestly why I gave as much preface as I did).

%% Simplified: In this case it would complain about "Unrecognized function or variable 'vals'.
function [odefunc] = funcGen(N)
odefunc = @(t,vals) [];
for i = 1:N
    n = 2*(i-1);
    recale = vals([n+1,n+2]) / norm(vals([n+1,n+2]));
    odefunc(t,vals) = @(t,vals) [odefunc(t,vals);vals(n+1);vals(n+2)];
end
end

[Michael Pilgrim]

%% Full Code
[odefunc, IC] = electricSet(10, 1e-12);
tspan = [0 1e-12];
[time, data] = ode45(odefunc,tspan,IC);



Local functions
function [odefunc, ICs] = electricSet(N, spacing)

pc = physConst();
c = pc.c;
q = [pc.q -pc.q];
m = [pc.mp pc.me];
elecprots = round(rand(1,N)+1,0); % Generates N 1's or 2's. 1's are protons and 2's are electrons.
q = q(elecprots);
m = m(elecprots);

ICs = zeros([1,6*N]);
for i = 1:N
    n = i-1;
    x = rand(1,1)*spacing;
    y = rand(1,1)*spacing;
    z = rand(1,1)*spacing;
    v = (.2*randn([1,1])+.4)*c;
    direction = rand(1,3)-.5;
    direction = direction / norm(direction);
    vx = v*direction(1);
    vy = v*direction(2);
    vz = v*direction(3);
    ICs(6*n+1:6*n+6) = [x y z vx vy vz];
end

[Ex, Ey, Ez] = electricField(N,q);

odefunc = @(t,vals) 0;
for i = 1:N
    n = i-1;
    
    odefunc = @(t,vals) [odefunc(t,vals);... % Repeating pattern of [x y z vx vy vz]
        vals(6*n+4);...
        vals(6*n+5);...
        vals(6*n+6);...
        Ex(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i);...
        Ey(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i);...
        Ez(vals(6*n+1),vals(6*n+2),6*n+3,vals)*q(i)/m(i)];
end

end

function [Ex, Ey, Ez] = electricField(N,q)
pc = physConst;
k = pc.kC;

Ex = @(x,y,z,vals) 0;
Ey = @(x,y,z,vals) 0;
Ez = @(x,y,z,vals) 0;

for i = 1:N
    n = i-1;
    r = norm([x-vals(3*n+1), y-vals(3*n+2), z-vals(3*n+3)]);
    direction = [x-vals(3*n+1), y-vals(3*n+2), z-vals(3*n+3)] / r;
    magnitude = k*q(i)/(r^2);
    Ex = @(x,y,z,vals) Ex(x,y,z,vals) + magnitude*direction(1);
    Ey = @(x,y,z,vals) Ey(x,y,z,vals) + magnitude*direction(2);
    Ez = @(x,y,z,vals) Ez(x,y,z,vals) + magnitude*direction(3);
end
end

[Michael Pilgrim]

Unrecognized function or variable 'x'.

Error in untitled2>electricField (line 58)

r = norm([x-vals(3*n+1), y-vals(3*n+2), z-vals(3*n+3)]);

Error in untitled2>electricSet (line 31)

[Ex, Ey, Ez] = electricField(N,q);

[Ameer Hamza]

In the simplified example, can you explain what are you trying to do with variable 'recale'? And what do you expect the value of 'val'? Note that 'val' is only defined in the scope of the anonymous function, i.e., only on that line. It does not have any meaning other than that line of code.

Also, how is physConst() defined?

[Walter Roberson]

    odefunc(t,vals) = @(t,vals) [odefunc(t,vals);vals(n+1);vals(n+2)];

You would have to use

   odefunc = @(t,vals) [odefunc(t,vals);vals(n+1);vals(n+2)];

At each step you would be shadowing the previous odefunc, "capturing" it inside the anonymous function, so that you would have a chain of anonymous functions that called previous versions of the anonymous function, and it's turtles all the way down to the original [] version. This is permitted, but it is difficult to debug, and it is inefficient. It does not build up a single anonymous function that creates all of the elements in a single step: it would be function calling function calling function... etc., until you got to the [] .