How to efficiently access surrounding matrix elements inside an ODE45 function.

I am working on a 3D simulation of magnetic moments which creates a problem similar to the Ising model: a problem of spins on a lattice. In order to find the lowest energy state, I am solving an ODE of type d theta / dt = f(theta, theta_i, t) which is a function of 8 surrounding spins theta_i (due to the shape of the unit cell). Currently the function I am using is

function [dthetas] = ODEsystemsMoments3DRealSymmetry(thetas,moments,K,A,C,n,m,o)
dthetas=zeros(o+2,m+2,n+2);

mythetas = reshape(thetas, o+2, m+2, n+2);
for i=2:o+1
    for j = 2:m+1
        for k = 2:n+1
            anisotropy = +K*cos(mythetas(i,j,k))*sin(mythetas(i,j,k));
            %disp(num2str(k))
            above1 = -A(i-1,j-1,k)*moments(i,j,k)*moments(i+1,j,k+1)*sin(mythetas(i,j,k) - mythetas(i+1,j,k+1));
            above2 = -A(i-1,j-1,k)*moments(i,j,k)*moments(i-1,j,k+1)*sin(mythetas(i,j,k) - mythetas(i-1,j,k+1));
            above3 = -A(i-1,j-1,k)*moments(i,j,k)*moments(i,j+1,k+1)*sin(mythetas(i,j,k) - mythetas(i,j+1,k+1));
            above4 = -A(i-1,j-1,k)*moments(i,j,k)*moments(i,j-1,k+1)*sin(mythetas(i,j,k) - mythetas(i,j-1,k+1));
            below1 = -A(i-1,j-1,k-1)*moments(i,j,k)*moments(i+1,j,k-1)*sin(mythetas(i,j,k) - mythetas(i+1,j,k-1));
            below2 = -A(i-1,j-1,k-1)*moments(i,j,k)*moments(i-1,j,k-1)*sin(mythetas(i,j,k) - mythetas(i-1,j,k-1));
            below3 = -A(i-1,j-1,k-1)*moments(i,j,k)*moments(i,j+1,k-1)*sin(mythetas(i,j,k) - mythetas(i,j+1,k-1));
            below4 = -A(i-1,j-1,k-1)*moments(i,j,k)*moments(i,j-1,k-1)*sin(mythetas(i,j,k) - mythetas(i,j-1,k-1));
            dthetas(i,j,k) = -(anisotropy + above1 + above2 + above3 + above4 + below1 +below2 + below3 + below4);
        end
    end
end
dthetas = reshape(dthetas, (n+2)*(m+2)*(o+2), 1);

and this works but is very slow (for 100x1000x64, it takes more than 3 days on an 8th gen i5).

I wondered if there is a way to improve the speed with gpuArrays or an efficient way to vectorise it but having to make arrays so regularly (because the angles update after each iteration) slows the function down drastically.

Now I am wondering if there is a nearest neighbour search or a similar tool that could efficiently extract the required moment and theta values.

Thanks very much for any insght.