Attempt to vectorize a Runge-Kutta algorith

Good morning everyone !

I'm currently playing with a code which consists of a scores of Matlab functions, and one of which relies on a Runge-Kutta algorith. I'd like to improve the length of time required to implement it, namely by vectorizing the already-existing iterative one.

I have been carrying out several attempts for more than a week, and unfortunately, the final result appears to be a dismal failure. If you catch a glimpse at the code written just after my comment, you could understand that my major issue remains the vectorization of 'statp'. Indeed, 'statp' takes part in the outter-loop, but all its parameters are implemented in the inner-loop.

Therefore, could you help me solving this problem by vectorizing this code ?

Here is the sketch of code : (tab, tbh_beta_abs, beta_abs, p0, phi, alfa dtanlam and curv have already been defined in other functions. pchip enables to carry out some interpolations)

gamma=1.4;
am=tab(1,:);
bm=tab(2,:);
len=sqrt((am-am(1)).^2+(bm-bm(1)).^2);
tub=len/(len(end));

statp(1)=10^5;

beta_abs_tub=pchip(tbh_beat_abs,beta_abs,tub);
f(1)=0;

for k=1:(length(tub)-1)

     dh=tub(k+1)-tub(k);
     dz=am(k+1)_am(k);
     dr=bm(k+1)-bm(k);
     dlen=len(k+1)-len(k);

     for j=1:4

         Z=flor(j/2);
         hh=tub(k)+.5*dh*Z;
         zz=am(k)+.5*dz*Z;
         rr=bm(k)+.5*dr*Z;

         statpp=statp(k)+.5*f(j)*Z;
         p0_hhm=pchip(tub,p0,hh);

         squaremak=((statpp/p0_hhm)^(-(gamma-1)/gamma)-1)*2/(gamma-1);

         newphi=pchip(tub,phi,hh);
         cosalfa=cos(alfa);
         sinalfa=sin(alfa);

         coslam=cos(newphi-alfa);
         sinlam=sin(newphi-alfa);
         tanlam=tan(newphi-alfa);

         dtanlam_int=pchip(tub,dtanlam,hh);

         curvature=pchip(tub,curv,hh);
         beta_abs_hh=pchip(tbh_beta_abs,beta_abs,hh);
         squaremaktan=(sqrt(squaremak)*sin(beta_abs_hh))^2;
         squaremaketa=(sqrt(squaremak)*cos(newphi-alfa)*cos(beta_abs_hh))^2;

         f(j+1)=gamma*statpp*dlen/(1-squaremaketa/coslam^2)*    ((1-squaremaketa)*cosalfa*squaremaktan/rr+...+curvature/(coslam^3)); %I have shortened the formula

         end

statp((k+1)=statp(k)+1/6*(f(2)+2*(f(3)+f(4))+f(5));

end

Thank you for your help !

[Oleg Komarov]

Please provide some example inputs so that we can run and profile it for you.

[Geoff]

I don't remember asking this question! ;-)