Speeding up loop-rich Matlab function to calculate temperature distribution I would like to speed up this function as much as possible in Matlab.
This is part of a bigger simulation project, and as it is one of the most called functions within the simulation, this is crucial.

For now, I tried generating a MEX file, but the speed was not better.

Vectorizing seems difficult (but would be beneficial due to the nested loops), given the non-linear operations.

What I'm trying to compute with this is the degree of temperature mixing between water layers, given by these equations:

$$\begin{align}
T_j =&amp;\ r\ T_{j+1} + (1-r)\ T_j\quad&amp;\textrm{for } j \in [i,\ i+\Delta i]\\
T_{j+1} =&amp;\ (1-r)\ T_{j+1} + r\ T_j
\end{align}$$

and \$n\$ is the total number of water layers. This must be calculated for all \$j\$, this is for all the layers, at all timesteps.

function y = mixing(T,dis,rr,n)  
%% =================================================================== 
%  input: temperature of cells array T, distance array dis, number of cells
%         n, mixing ratio r
%
%  output:  new temperature array
%
%  purpose: calculates the temperature array of next timestep
%  ===================================================================

for j = 1:n
    i = 1;
    r = rr;
    while i &lt; dis(j)+1 &amp;&amp; j+i &lt;= n
        if (dis(j) &lt; i) 
            r = r*(dis(j)-floor(dis(j))); 
        end
        d = T(j+i-1);
        T(j+i-1) = r*T(j+i) + (1-r)*T(j+i-1);
        T(j+i) = r*d + (1-r)*T(j+i);
        i = i + 1;
    end
end

y = T;
end


Inputs: T is a 10 × 1 double, dis is a 10 × 1 double, rr is a 1 × 1 double, and n is a 1 × 1 integer value.
First, using vectors would probably speed it up quite much, and would be worth the effort.

However, here's a small tip, that could improve the performance a little bit:

In many of the CPU-hungry iterative computations, you could replace the computations that occur several times by a variable. 

In the code, there are several "j+i" terms. They can simply be replaced by a variable named, "jpi". So, the number of additions in each value of j would decrease by 5 (not 6, as you still need to do an addition, once per each iteration of j).

dis(j) is also frequently used, and it might be useful to create another dummy variable, like dis_j = dis(j) once per each iteration of j. I don't know if matlab compiler does this itself to optimize performance. this one stops the cpu from looking at the value of dis(j) each time it is called. So it should help a little bit.