How to make a function that uses Runge-Kutta Method

Here is my code to find the velocity and position of particles in 3 dimensions using Rung-Kutta 4th order method as a time marching:

I just mentined a part of code here,its a little long.

Uf=Poiseuille flow(fluid) in pipe and wp,up,vp are particles velocity in 3 dimensions and FW=dwp/dt,FU=dup/dt,FV=dvp/dt.

k1-k4 and m1-m4 and l1-l4 are slopes from Runge-Kutta 4th oreder method to obtain velocity and position of particles.

MY QUESTION is how can I read a Runge-Kutta 4th as an function to avoid repeating these k,m and l's and just call the that function.

while ( t(j)<Time )
       FW=@(wp,Uf) (1/(rop+0.5*ro))*((rop-ro)*g+(18*miu/Dp^2)*(Uf-wp));  % wp is z's particle velocity
       FU=@(up,Uf) (1/(rop+0.5*ro))*((rop-ro)*0+(18*miu/Dp^2)*(0 -up));  % up is x's particle velocity
       FV=@(vp,Uf) (1/(rop+0.5*ro))*((rop-ro)*0+(18*miu/Dp^2)*(0 -vp));  % vp is y's particle velocity
         k1 = FW(wp(:,:,:,j),Uf)*Dt;        m1 = FU(up(:,:,:,j),Uf)*Dt;        l1 = FV(vp(:,:,:,j),Uf)*Dt;
         k2 = FW(wp(:,:,:,j)+0.5*k1,Uf)*Dt; m2 = FU(up(:,:,:,j)+0.5*m1,Uf)*Dt; l2 = FV(vp(:,:,:,j)+0.5*l1,Uf)*Dt;
         k3 = FW(wp(:,:,:,j)+0.5*k2,Uf)*Dt; m3 = FU(up(:,:,:,j)+0.5*m2,Uf)*Dt; l3 = FV(vp(:,:,:,j)+0.5*l2,Uf)*Dt;
         k4 = FW(wp(:,:,:,j)+k3,Uf)*Dt;     m4 = FU(up(:,:,:,j)+m3,Uf)*Dt;     l4 = FV(vp(:,:,:,j)+l3,Uf)*Dt;
         
         wp(:,:,:,j+1) = wp(:,:,:,j) + (1/6)*(k1+2*k2+2*k3+k4);   % particle velocity along the pipe.
         up(:,:,:,j+1) = up(:,:,:,j) + (1/6)*(m1+2*m2+2*m3+m4);
         vp(:,:,:,j+1) = vp(:,:,:,j) + (1/6)*(l1+2*l2+2*l3+l4);
         
         Zp(:,:,:,j+1) = Zp(:,:,:,j)+Dt*wp(:,:,:,j);              % particle displacement along the pipe.
         Xp(:,:,:,j+1) = Xp(:,:,:,j)+Dt*up(:,:,:,j);
         Yp(:,:,:,j+1) = Yp(:,:,:,j)+Dt*vp(:,:,:,j);
         
         j=j+1;
         t(j)=Dt*j;
end

Answer 1

[darova]

Try something like that

wp(j+1) = rk4(FW,t,wp(j),h);
function u1 = rk4(f,t,u0,h)
    k1 = f(t,u0);
    k2 = f(t+h/2,u0+k1*h/2);
    k3 = f(t+h/2,u0+k2*h/2);
    k4 = f(t+h,u0+h*k3);
    u1 = u0 + 1/6*h*(k1 + 2*k2 + 2*k3 + k4);

[H-M]

Thank you Darova for your answer

[darova]

my pleasure

[H-M]

Darova,Do I need to add for loop for iteration or not?

[darova]

I think so

[H-M]

Darova would you please give me an idea to write the for loop, it starts from i=1 to what?

thank you

for i=1:(what?)

[darova]

i think something like that

for i = 1:length(t)
    wp(j+1) = rk4(FW,t(i),wp(i),Dt);
end