How can I make my code run faster?

a = 1; 
b = 1;
rho = 1;
D = 3;
h =1;

%set matrix phi
Nx=150;
Ny=100;
phi = rand(Nx,Ny);
%assign random -1 and 1 for phi matrix
for i=1:Nx
    for j=1:Ny
        if phi(i,j)<=0.5
            phi(i,j)=-1;
        else
            phi(i,j)=1;
        end
    end
end
Problem = input ('problem number');
switch Problem
   case 1
  %Five Point Stencil
dt = 1e-3; 
for  k = 1:dt:20    
 % apply Laplacian to matrix phi   
 v = Laplacian_2D (phi,h,5);
 % Find MAtrix A in equation 4  
 A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
 % Calculate Laplacian of Matrix A for equation 4
 %left hand of equation 4
 E4 = D * Laplacian_2D (A,h,5);
 %First Order Runge Kutta with Equation 4
 c1 = dt*E4;
 phi = phi + c1;
%draw plot
 imagesc (phi);
 colorbar;
 caxis ([-1,1])
 title ('Five Point Stencil: First Order Runge Kutta');
 drawnow;
 end
    case 2
    %Nine Point Stencil
 dt = 1e-4;
 for k = 1:dt:(20/dt)
  counter = 0;
  % apply Laplacian to matrix phi   
 v = Laplacian_2D (phi,h,9);
 % Find MAtrix A in equation 4  
 A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
 % Calculate Laplacian of Matrix A for equation 4
 v2 = Laplacian_2D (A,h,9);
 %left hand of equation 4
 E4 = D * v2;
 %First Order Runge Kutta with Equation 4
 c1 = dt*E4;
 phi = phi + c1;
%draw plot
 imagesc (phi);
 colorbar;
 caxis ([-1,1])
 title ('Nine Point Stencil: First Order Runge Kutta');
 drawnow;
 
 end
 
    case 3
    %Second Order Runge Kutta with Five Point Stencil
   dt = 5e-3;
for  k = 1:dt:(20/dt)
 
 % apply Laplacian to matrix phi   
 v = Laplacian_2D (phi,h,5);
 % Find MAtrix A in equation 4  
 A = ((b^4*phi.^3 - a*b^2*phi)-rho*v);
 % Calculate Laplacian of Matrix A for equation 4
 v2 = Laplacian_2D (A,h,5);
 %left hand of equation 4
 E4 = D * v2;
 %run again with c1
 c1 = dt*E4;
 v3 = Laplacian_2D (phi+0.5*c1,h,5);
 R = ((b^4*(phi+0.5*c1).^3 - a*b^2*(phi+0.5*c1))-rho*v3);
 v4 = Laplacian_2D (R,h,5);
 E4_2 = D * v4;
 %Second Order Runge Kutta with Equation 4
 c2 = dt*E4_2;
 phi = phi + c2;
 %draw plot
 imagesc (phi);
 colorbar;
 caxis ([-1,1])
 title ('Five Point Stencil: Second Order Runge Kutta');
 drawnow;
end
end

function v = Laplacian_2D(u,h,stencil)
Nx = 150; 
Ny = 100;
v = zeros (Nx,Ny); 
if stencil == 5 %insert for the five point stencil
%periodic boundary conditions
    for i = 1:Nx
            N = i-1;
            if N==0
               N=Nx;
            end            
            S = i + 1;
            if S== Nx+1
                S=1;
            end   
    for j = 1:Ny       
            E = j+1;
            if E == Ny+1
                E = 1;
            end            
            W = j-1;
            if W == 0
                W=Ny;
            end 
             v (i,j) = (1/h^2)*(u(N,j) + u(S,j) + u(i,W) + u(i,E) - 4*u(i,j));
    end
    end          
elseif stencil == 9 %insert for the nnine point stencil
   %periodic boundary conditions 
    for i = 1:Nx
         N = i-1;
            if N==0
               N =Nx;
            end            
            S = i + 1;
            if S== Nx+1
                S = 1;
            end
    for j = 1:Ny                     
            E = j+1;
            if E == Ny+1
                E = 1;
            end
            W = j-1;
            if W == 0
                W=Ny;
            end 
             v (i,j) = (1/(6*h^2))*(u(N,W)+ 4*u(N,j) + u(N,E) + 4*u(i,W) + 4*u(i,E) + u(S,W) + 4*u(S,j) + u(S,E) - 20*u(i,j));
    end
    end
else 
    error ('stencil must be 5 or 9');
end
end

[per isakson]

Am I right that the purpose of the script is to create a picture of how phi develops ? With a faster code this will be even more difficult to observe.

[Rik]

Have you run the profiler to find the bottleneck?

[Adam]

First step in working out how to make code faster is to understand what is slow, as Rik says.

Answer 1

[darova]

Shorter version of Laplacian_2D function. Is there any connection to gradient and laplacian built-n functions?

un = circshift(u, [-1 0]);
us = circshift(u, [+1 0]);
uw = circshift(u, [0 -1]);
ue = circshift(u, [0 +1]);

unw = circshift(u, [-1 -1]);
une = circshift(u, [-1 +1]);
usw = circshift(u, [+1 -1]);
use = circshift(u, [+1 +1]);

if stencil == 5
    v = 1/h^2*(un + us + uw + ue - 4*u);
elseif stencil == 9 %insert for the nnine point stencil
    v = 1/(6*h^2)*(4*(un+uw+ue+us) + unw+une+usw+use - 20*u);
else 
    error ('stencil must be 5 or 9');
end

This part can be shorter

% phi = rand(Nx,Ny);
% %assign random -1 and 1 for phi matrix
% for i=1:Nx
%     for j=1:Ny
%         if phi(i,j)<=0.5
%             phi(i,j)=-1;
%         else
%             phi(i,j)=1;
%         end
%     end
% end
phi = randi([0 1],Nx,Ny)*2-1;