And some cleanup:
...
r = (0:N_i) * dr;
c = zeros(N_i, N_t - 1);
c(:, 1) = c_inf;
c(r < a_0, 1) = c_s;
...
dr3 = dr ^ 3;
...
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
end
Best Answer