This code is taking a long time to finish. I tried with vectorizing the inner loops n=2:ny-1, m=2:nx-1, but it is still slow. Should I also vectorize the time loop?

nx=320;
ny=320;
nt=4000;
predictor = zeros(nt,nx,ny);
corrector = zeros(nt,nx,ny);
c4 = zeros(nt,nx,ny);  
%there are some initial conditions which I did not mention here
for t=1:nt           % As the scalar starts to get updated from t=2, we continue the loop until t=(nt+1) 
    for n=2:ny-1
        for m=2:nx-1
            predictor(t,m,n)=c4(t-1,m,n)-((u*dt)/dx)*(c4(t-1,m+1,n)-c4(t-1,m,n))-((v*dt)/dy)*(c4(t-1,m,n+1)-c4(t-1,m,n));
            corrector(t,m,n)=0.5*(c4(t-1,m,n)+predictor(t,m,n)-((u*dt)/dx)*(predictor(t,m,n)-predictor(t,m-1,n))-((v*dt)/dy)*(predictor(t,m,n)-predictor(t,m,n-1)));
        end
    end
    c4=corrector;
    
    % Visualizing the evolution of the scalar with time (at each time)
    set(0,'CurrentFigure',f4);
    [X,Y]=meshgrid(0:dx:(nx-1)*dx,0:dy:(ny-1)*dy);
    surf(X,Y,squeeze(c4(t,1:nx,1:ny)));
    shading flat;
    colorbar;
    drawnow;
end

Answer 1

[Matt J]

Yes, there is no reason to be doing any looping when constructing the variables predictor and corrector. They can be computed in one line each with vectorized arithmetic operations. Once you've computed these, you can run a loop over t to make your surface plots, but that's the only thing in your code that requires a loop.

Also, you should probably not be plotting 4000 surfaces. That will take a long time to render no matter how you do it. Maybe plot a subset of the full set.

Also, you should also definitely not be doing this inside the loop,

[X,Y]=meshgrid(0:dx:(nx-1)*dx,0:dy:(ny-1)*dy);

You are simply repeating the same (expensive) operation 4000 times. You can do this once prior to the loop or even not at all, since surf() does not require the X and Y inputs to be meshgrid matrices. They can be simple vectors.

    [X,Y]=deal(0:dx:(nx-1)*dx,0:dy:(ny-1)*dy);
    
    for t=1:40:nt
     surf(X,Y,squeeze(c4(t,:,:))); drawnow
    end