MATLAB: Caught in loop

loops

My code is as follows:
alpha = 0.5;
beta = 0.1;
a_0 = 1;
dr = 0.1;
dt = input('Value of dt:');
h = 0.001;
c_inf = 0;
c_a = 1;
c_s = 2;
mesh_points = 3;
rho = dt/(dr)^2;
R= 30;
N_i = R/dr;
t_diss = 5.8;
N_t = t_diss/dt;
if rho >= 0.5;
    disp('Invalid dt value');
    return;
end
for i=1:N_i;
    r(i)=(i-1)*dr;
    if r(i)< a_0
        c(i,1)= c_s;
    else 
        c(i,1)= c_inf;
    end
end
a(1) = a_0;
for j = 1:N_t -1;
    b_1 = floor (a(j)/dr) + 2;
    b_2 = b_1 + 1;
    c_ah = (a(j)+h-b_1)*(a(j)+h-b_2)*(c_a)/((a(j)-b_1)*(a(j)-b_2))+h*(a(j)+h-b_2)*c(b_1,j)/((b_1-a(j))*(b_1-b_2)) + h*(a(j)+h-b_1)*c(b_2,j)/((b_2-a(j))*(b_2-b_1));
    dc_dr = ((c_ah)-c_a)/h;
    a(j+1) = a(j) +dt*beta*dc_dr;
     for i=2:N_i - 1;
         if r(i) > a(j+1);
             d = rho - (rho/(i-1)) - (a(j)^2)*(alpha-1)*(a(j+1)-a(j))/(2*(i-1)^2*(dr)^3);
             e = 1-2*rho;
             f = rho + (rho/(i-1)) + (a(j)^2)*(alpha-1)*(a(j+1)-a(j))/(2*(i-1)^2*(dr)^3);
             c(i, j+1) = d*c(i-1,j) + e*c(i,j) + f*c(i+1, j);
         elseif r(i) == a(j+1);
                 c(i,j+1) = c_a;
             else c(i,j+1) = c_s;
         end
             disp(c(i,j));
         end
         c(1, j+1) = c(2, j+1);
  end
It basically calculates concentration near a dissolving sphere but when it gets to disp(c(i,j)) I get caught in loop of displaying numbers that change periodically. Any ideas?

Best Answer

And some cleanup:
...
% Vectorization of the "for i=1:N_i;" loop:
r = (0:N_i) * dr;
c = zeros(N_i, N_t - 1);  % Pre-allocate!
c(:, 1)       = c_inf;
c(r < a_0, 1) = c_s;
...
dr3 = dr ^ 3;
...
  % Reduce repeated calculations inside the "for j = 1:N_t -1" loop:
  ajp1   = a(j) +dt*beta*dc_dr;
  a(j+1) = ajp1;
  tmp1   = (a(j)^2) * (alpha-1) * (ajp1-a(j)) / (2 * dr3);
  for i = 2:N_i - 1
    if r(i) > ajp1
      tmp2 = rho / (i-1) + tmp1 / ((i-1)^2);
      d = rho - tmp2;
      e = 1 - 2 * rho;
      f = rho + tmp2;
      c(i, j+1) = d*c(i-1,j) + e*c(i,j) + f*c(i+1, j);
    elseif r(i) == ajp1
      c(i, j+1) = c_a;
    else
      c(i, j+1) = c_s;
    end
    % disp(c(i,j));
  end
Related Question