Reducing run time of a numeric calculation using a mex file in Matlab

I wrote a Matlab code that involves doing a numeric calculation (relaxation), but it is quite slow. I learned of the possibility of using a mex file to run a C code and integrate it into Matlab, so I was thinking of doing the numerical calculation (which is relatively simple but involves loops and takes time) in C, and the rest (before and after) in Matlab.

The part of my Matlab code where the calculation is done:

% evolution of the potentials %
% note: for index directions with periodic boundary conditions: index=mod(index-1,L)+1.
%     for index=index+1 it is mod(index,L)+1 , 
% and for index=index-1 it is mod(index-2,L)+1
for i_t=1:max_relaxation_iterations
  for q=1:length(i_eff_V_bounded) % this is set instead of running i=2:(L-1), j=1:L , k=1:L and ending up going over sites that are 0 in our effective system %
    i=i_eff_V_bounded(q);
    j=j_eff_V_bounded(q);
    k=k_eff_V_bounded(q);
    V0=V(i,j,k);    
    V1=( V(i+1,j,k)+V(i-1,j,k)+V(i,mod(j,L)+1,k)+V(i,mod(j-2,L)+1,k)+V(i,j,mod(k,L)+1)+V(i,j,mod(k-2,L)+1) )/( system(i+1,j,k)+system(i-1,j,k)+system(i,mod(j,L)+1,k)+system(i,mod(j-2,L)+1,k)+system(i,j,mod(k,L)+1)+system(i,j,mod(k-2,L)+1) ); % evolving the potential as the average of its occupied neighbors %
    V(i,j,k)=V0+(V1-V0)*over_relaxation_factor; % evolving the potentials in time with the over relaxation factor %
    delta_V_rms(i_t)=delta_V_rms(i_t)+(V1-V0)^2; % for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg %
    delta_V_abs(i_t)=delta_V_abs(i_t)+abs(V1-V0); % for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg %                
    delta_V_max(i_t)=max(abs(V1-V0),delta_V_max(i_t)); % for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg % 
  end  
end

So in C it should be something like:

#include <stdio.h>

int mod(int x,int N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */

{

return (x%N+N)%N;

}

double d_abs(double x) /* a function for the absolute value operator */

{

if x<0

{

return -x;

}

else

{

return x;

}

}

double max(double x,double y) /* a function for the max operator */

{

if x>y

{

return x;

}

else

{

return y;

}

}

/* evolution of the potentials /*

/* note : periodic boundary conditions for the j,k directions */

void potentials_evolution(int max_relax_iters,int N_eff_occ_sites,int i_eff_V_bounded[],int j_eff_V_bounded[],int k_eff_V_bounded[],int system[][][],over_relax_fact,double V[][][],double delta_V_rms[],double delta_V_abs[],double delta_V_max[])

{

int i_t,q,i,j,k;

double V0,V1;

for(i_t=0;i_t<max_relax_iters;i_t++)

{

for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */

{

i=i_eff_V_bounded[q];

j=j_eff_V_bounded[q];

k=k_eff_V_bounded[q];

V0=V[i][j][k];

V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */

V[i][j][k]=V0+(V1-V0)*over_relax_fact; % evolving the potentials in time with the over relaxation factor %

delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */

delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */

delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */

}

}

}

And so in Matlab I will replace the part of my Matlab code shown above with something like:

potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);

How do I implement this? I tried looking for a simple way to do it but I couldn't figure how to properly do it.

Note 1: This numeric calculation is done not just once but many times for different systems that are generated randomly (there is a `for` loop going over the different systems).

Note 2: I've never used mex files and my C is quite rusty.