Question on implementing Danckwert's boundary condition

I have the following code, which solves for the change in concentration of a species along the length of a channel of length L . The equation is ,

$ \frac{\partial C}{\partial t} = -v\frac{\partial C}{\partial z} + D\frac{\partial^2C}{\partial z^2} - kC $

The spatial direction is discretized and the resulting form solved using ode solver.

I have used constant values of concentration at the right and left boundary of the cylindrical channel.

    function Model()
    Nodes = 10;
    v = 10;D =50;k = 10;
    %C0 contains the concentration of the species at time 0 at all the nodes 
    C0(1:10) = [2000;1996;1980;1973;1965;1950;1947;1938;1878;1800]';
    C(1:8,1) = C0(2:9);

    tspan  = 1:10;
    [t C]  = ode15s(@(t,v) diff(t,v), tspan ,C)
    function dCdot=diff(t,C)
    Convection = zeros(8,1);
    ConvectionB = zeros(8,1);
    ConvectionB(8,1) = C0(10);%Value at the right boundary(z=L) goes here
    Matrix = -1*diag(ones(1,8));
    for i=1:7
        Matrix(i,i+1) = 1;
    end
    Convection  = v*(Matrix*C + ConvectionB);
    %The second term on the LHS of the first equation is discretized using 
    %central difference and a tridiagonal matrix is formed
    Diffusion = zeros(8,1);
    DiffusionB  = zeros(8,1);
    DiffusionB(1,1) = C0(1);%Value at the left boundary(z=0) goes here
    DiffucionB(8,1) = C0(10);%Value at the right boundary(z=L) goes here
    Diffusion  =  D*(full(gallery('tridiag',8,1,-2,1))*C + DiffusionB);

    Reaction = -k*C;

    dCdot = Convection + Diffusion + Reaction; 
    end

    end

I would like to implement Danckwert’s boundary condition given as follows,

At z = 0,

C (z =0) = C^{feed} + \frac{D}{v}\frac{\partial C}{\partial z}

At z = L,

\frac{\partial C}{\partial z} = 0

Could someone advise me on how the above boundary condition can be implemented in the code?

[Ngo Ich Son]

You only can apply the Danckwert's BC using the boundary value problem solver (bvp4c or bvp5c).