How to code Flux Boundary Conditions for a PDE using ODE15s and method of lines

I have a PDE in spherical coordinates that I've discretized spatially using the method of lines for solving with ODE15s. I am having trouble figuring out how to represent the flux boundary conditions:

1.) At r=0, ∂c/∂r = 0

2.) at r=R, De ∂c/∂r = kf (C − c)

My code looks like this: Thanks for your help!

function [dy] = Kinetics(t,y)
kf=1.8*10^-3; %Mass transfer coefficient (cm/s)
De=8*10^-8;
R=20*10^-4;
ka=2*10^-4;
Keq=72;
Vm=40;
V=40;
n=100;
qm=45;
Ep=0.6;
cp=y(1:n);
Qp=y(n+1:2*n);
Cp=y(2*n+1);
dCp=zeros(1,1);
dcp=zeros(n,1);
dQp=zeros(n,1);
Vl_stage=(1/n);

%%At r=0 Boundary Conditions
dcp(1)=(cp(1)-cp(2))/Vl_stage; %This is where I tried to implement Boundary Condition at r=0
dQp(1)=((ka)*(((qm-Qp(1))*cp(1))-(Qp(1)/Keq)));


for i=2:n-1
    Ri=(Vl_stage*(i)/R);
    Ri1=(Vl_stage*(i+1)/R);
    dcp(i)=(De/Ep)*(1/Ri1^2)*(((Ri1^2*(cp(i+1)-cp(i)))-(Ri^2*(cp(i)-cp(i-1))))/(Vl_stage^2))-(dQp(i)/Ep);
    dQp(i)=((ka)*(((qm-Qp(i))*cp(i))-(Qp(i)/Keq)));
end

%%At r=R Boundar Conditions
for i=n  
    dcp(n)=(kf/De)*(Cp(1)-cp(n))/Vl_stage; %This is where I tried to implement Boundary Condition at r=R
    dQp(n)=((ka/Ep)*(((qm-Qp(n))*cp(n))-(Qp(n)/Keq)));
    dCp(1)=-(3*Vm/(R*V))*kf*(Cp(1)-cp(n));
end


dy=[dcp;dQp;dCp];

end

And I am calling with [t,y] = ode15s(@Kinetics,[0,1],[horzcat(zeros(1,200),10)]);

If it helps, here is the original system of equations:

Answer 1

[Torsten]

One of many possibilities:

dc/dr = 0

c2 - c1 = 0

->

dc1/dt = dc2/dt

De*dc/dr = kf*(C-c)

De*(cn-c(n-1))/Vl_stage = kf*(C-cn)

cn = (kf*C + De*c(n-1)/Vl_stage)/(kf+De/Vl_stage)

->

dcn/dt = (kf*(dC/dt) + De*(dc(n-1)/dt)/Vl_stage)/(kf+De/Vl_stage)