Rungee Kutta Example

The following code is just a trial feedback is much appriciated to check the code, a note to the user that this a computationaly heavy code to run so you might need to use a more coarse grid.

The equations will be written here:

clc
clear
M=32;
N=32;
L=32;
LX=1;
LY=1;
LZ=1;
dx=LX/M;
dy=LY/N
dz=LZ/L
Tt=1;
t0=0;
T=50;
dt=(Tt-t0)/T;
u0(1:M,1:N,1:L,1)=0;
v0(1:M,1:N,1:L,1)=0;
w0(1:M,1:N,1:L,1)=0;
for t=1:T;
for i=1:M;
    for j=1:N;
        for k=1:L;
            u(i,j,k,1)=u0(i,j,k,1);
            v(i,j,k,1)=v0(i,j,k,1);
            w(i,j,k,1)=w0(i,j,k,1);
            uRHS(i,j,k,t)=randn(1,1);
            vRHS(i,j,k,t)=randn(1,1);
            wRHS(i,j,k,t)=randn(1,1);
            k1x(i,j,k,t)=dt*uRHS(i,j,k,t);
            k2x(i,j,k,t)=dt*(0.5*k1x(i,j,k,t)+dt*uRHS(i,j,k,t));
            k3x(i,j,k,t)=dt*(0.5*k2x(i,j,k,t)+dt*uRHS(i,j,k,t));
            k4x(i,j,k,t)=dt*(k3x(i,j,k,t)+dt*uRHS(i,j,k,t));
            u(i,j,k,t+1)=u(i,j,k,t)+(1/6)*k1x(i,j,k,t)+(1/3)*(k2x(i,j,k,t)+k3x(i,j,k,t))+(1/6)*k4x(i,j,k,t);
            k1y(i,j,k,t)=dt*vRHS(i,j,k,t);
            k2y(i,j,k,t)=dt*(0.5*k1y(i,j,k,t)+dt*vRHS(i,j,k,t));
            k3y(i,j,k,t)=dt*(0.5*k2y(i,j,k,t)+dt*vRHS(i,j,k,t));
            k4y(i,j,k,t)=dt*(k3y(i,j,k,t)+dt*vRHS(i,j,k,t));
            v(i,j,k,t+1)=v(i,j,k,t)+(1/6)*k1x(i,j,k,t)+(1/3)*(k2x(i,j,k,t)+k3x(i,j,k,t))+(1/6)*k4x(i,j,k,t);
            k1z(i,j,k,t)=dt*wRHS(i,j,k,t);
            k2z(i,j,k,t)=dt*(0.5*k1z(i,j,k,t)+dt*wRHS(i,j,k,t));
            k3z(i,j,k,t)=dt*(0.5*k2z(i,j,k,t)+dt*wRHS(i,j,k,t));
            k4z(i,j,k,t)=dt*(k3z(i,j,k,t)+dt*wRHS(i,j,k,t));
            w(i,j,k,t+1)=w(i,j,k,t)+(1/6)*k1z(i,j,k,t)+(1/3)*(k2z(i,j,k,t)+k3z(i,j,k,t))+(1/6)*k4z(i,j,k,t);
            v(i,j,k,t+1)=((u(i,j,k,t+1))^2+(v(i,j,k,t+1))^2+(w(i,j,k,t+1))^2)^0.5;
        end
    end
end
end
for t=1:T;
contour(v(:,:,3,t))
pause(0.2)
end

The code output: