MATLAB: Parallelized ODE45 solution with 3D spline

3d spline ode45 parfor

Hi all,
I have a set of data describing an electric potential distribution in space V(x,y,z) with 3 independent variables x, y, and z.
I want to use this function in solving a system of ODEs describing the 3D motion of a charged particle existing in this electric potential distribution.
I have no analytic expression for V(x,y,z) and the function is ill-fitted using polynomial fits, so I use a 3D spline to construct a piecewise polynomial to be used by the ODE45 solver.
Since the 3D spline introduces a large overhead time, I do it externally in a separate function then make the spline a global variable to be accessed by the ODE solver.
My problem is that now I want to solve the system of ODEs with different combinations of initial conditions and I want to parallelize the operation using parfor loops.
However, parfor loops are not allowing the use of global variables even though they are not changed within the loop itself.
Any suggestions on how I could parallelize the operation without forsaking the use of the global variables in reconstructing the function?
Here's the code mentioned:
global spx spy spz 
 sp=csapi({z_cord,y_cord,x_cord}, Data);
 spx=fnder(sp, [0,0,1]);
 spy=fnder(sp, [0,1,0]);
 spz=fnder(sp, [1,0,0]);
 
 
y0=[-12:0.01:-11];
for j=1:length(y0)
tstart=tic;
[t,sol]=ode45(@func_xy,[0 30e-6],[0;0;y0(j);0;0;0], odeset('Events', @(t,y) myEvent(t,y,tstart)));
end
function dxdt=func_xy(t,y)
global spx spy spz 
q=1.6e-19;
m=400*1.6e-27;
dxdt=[y(2);
    (-q/m)*fnval(spz,{y(1), y(3), y(5)}) ;
     y(4);
    (-q/m)*fnval(spy,{y(1), y(3), y(5)});
     y(6);
    (-q/m)*fnval(spx,{y(1),y(3), y(5)});];
end
function [value, isterminal, direction] = myEvent(t, y,tstart)
value      = double(toc(tstart)>200);
isterminal = 1;   % Stop the integration
direction  = 0;
end

Best Answer

Can you just pass the variables spx spy spz to func_xy as additional arguments?
q=1.6e-19;
m=400*1.6e-27;
odeOpts = odeset('Events', @(t,y) myEvent(t,y,tstart))
% pass additional arguments by using function handle:
odefunHandle = @(t,y)func_xy(t,y , q,m,spx,spy,spz)
[t,sol]=ode45(odefunHandle,[0 30e-6],[0;0;y0(j);0;0;0], odeOpts);
%%%
function dxdt=func_xy(t,y , q,m,spx,spy,spz)
% global spx spy spz 
% q=1.6e-19;
% m=400*1.6e-27;
dxdt=[y(2);
    (-q/m)*fnval(spz,{y(1), y(3), y(5)}) ;
     y(4);
    (-q/m)*fnval(spy,{y(1), y(3), y(5)});
     y(6);
    (-q/m)*fnval(spx,{y(1),y(3), y(5)});];
end
By the way are you sure that the main cost of the spline interp is the up-front cost of csapi/fnder, and not the fnval()s?
Related Question