Profiling your code it comes up that the most demanding routine is the Green function calculation, with 0.374s:
I have rewritten your function, eliminating the innermost for loop and vectorizing the code. I made explicit use of bsxfun even if, starting from Matlab2016b, it is possible to use implicit expansion. I prefer to clearly have under control what is happening in the code.
I also removed the last loop to update the final part of the greens function. The result is the following. I renamed the function as Greens_fn2 to avoid confusion (Remember to change the call in your u function).
function G = Greens_fn2(D,mu_n,r,t,r_bar)
G = G+bsxfun(@times,bsxfun(@times,(1+mu_n(i)^-2).*sin(mu_n(i)*r).*sin(mu_n(i)*r_bar),exp(-D*mu_n(i)^2*t(:))),(2*r_bar/r));
G = bsxfun(@plus,G,3*r_bar.^2);
A new profile, gives the new performance time:
that is 0.094s, with a speedup of about 4. The results are identical.
Possible additional measures to speedup the code:
- I am not sure but maybe you can also expand with respect to N_t, removing also the remaining loop.
- Use a mex-file for the calculation of the Green's function
Befor any other attempt, I would analyze how this code scales with respect to N_r and N_t (I think you are interested in more complex cases). Please let us know!