[exnihilo]

I am trying to integrate the system I am showing in the page 3 of the attached pdf but I am having some problems with Cw and the radius because they should reach a minimum and a maximum respectively but in my simulation they go towards +inf and -inf, and Cw is decreasing too fast. Cw0, K, Cm are parameters, I also attach the whole paper if something is not clear enough. (I am using matlab). I am using forward euler and finite difference method. It is probably better to integrate with a variable space step but I cant figure out how to do it, so I just made a big enough domain to allow the radius to increase. I have started recently cfd so be kind ahahah.

clc; clear all; close all
%% DATA
K = 1;
Diff = 4e-9;
Cw0 = 55e-3; %mol/m3
r0 = 1e-3;
L = 12e-3;
N = 500;
tf = 150;
Cm = 30e-3; %mol/m3
Cinf = 0;
%% preprocessing
h = L/N;
grid1D = linspace(0,L,N+1);
%% SOLUTION
index = find(grid1D >= r0);
dt = 0.005;
tsteps = tf/dt;
%initial sol
C = [Cw0.*ones(1,index(1)-1)./K,Cm*ones(1,index(end)-index(1))]';
%mem allocation for plot
Cplot = zeros(size(grid1D));
%loop
for t = 1:tsteps
    C0 = C;
    % bc 6: eq 6 integration
    if t ~= 1
        indexrd = find(grid1D >= rd);
        for counter_bc6 = index(1):indexrd(1) - 1
            int6(counter_bc6) = (-4*pi*C(counter_bc6)*grid1D(counter_bc6)^2 -4*pi*C(counter_bc6 + 1)*grid1D(counter_bc6+1)^2)*h/2;
        end
        Cw = sum(int6)/(4/3*pi*r0^3) + Cw0;
    else
        Cw = Cw0;
    end
    % eq 3
    C(index(1)) = Cw/K;
    % eq 5
    if t == 1
        rd = r0 + h;
    else
        d = -Cm + 8*Cm -8*C(indexrd(1)- 2) + C(indexrd(1)- 3);
        rd = -Diff*(d)/12/h/Cm*dt + rd;
    end
    if (mod(rd,h)~=0)
        rd = rd + (h - mod(rd,h));
    end
    %eq 4
    indexrd_new = find(grid1D>=rd);
    %expl
    for i = index(1):N-1
        C(i) = C0(i) + dt*(Diff/h^2*(C0(i+1) + C0(i-1) - 2*C0(i)) + 2*Diff/i/h*(C0(i+1)-C0(i-1))/h );
    end
    for i = indexrd_new(1) :length(C)
        C(i) = Cm;
    end
    for i = 1:index
        C(i) = Cw;
    end
    if (mod(t,100)==0)   % => Every 50 time steps
        for i=1:N-1
            Cplot(i) = 0.5*(C(i+1) + C(i));
        end
        plot(grid1D,Cplot);
        hold on
        plot(grid1D(indexrd_new(1))*ones(2,1),[0.028 0.052])
        %scatter(t*dt,Cw)
        hold off
        ylim([2.8e-2 5.2e-2])
        xlim([0 L/2])
        xlabel("length [m]")
        ylabel("Concentration [mol/m3]")
        message = sprintf('t=%d', ceil(t*dt));
        time = annotation('textbox',[0.15 0.8 0.1 0.1],'String',message,'EdgeColor','k','BackgroundColor','w');
        drawnow;
    end
end

.rtcContent { padding: 30px; }.lineNode {font-size: 12pt; font-family: "JetBrainsMono Nerd Font", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }

[hunt_mat]

The way to include Neumann boundary conditions into the diffusion is to use ghost cells at each point. You can then correctly fit in the ghost points that you get from the basic finite difference, that's how you do this normally. You have a moving boundary condition, so what you could do is fix the boundary with the transformation:

r is mapped to r/R(t), then this fixes the boundary. The resulting equations you get will be a coupled set of equations for the nodes and the moving boundary. You can then use Newton-Raphson to solve these equations if they're nonlinear, or a matrix method if they're linear.