function cwf_prod = ADE_model(qp,inj,T,cinj,NDiv,Lwell,Vp,D_ref,qf_ref)

cwf_prod = zeros(length(T),1);
inj = inj*5.615;
qp = qp*5.615;
Vp = 10^Vp;
q_ref = 0.5*(inj(1)+qp(1));
% D = 10^D;
qf_ref = 10^qf_ref;
dx = Lwell/NDiv;
C = zeros(NDiv,1);
Tn = 0;
for t = 1:(length(T))
    dt = T(t) - Tn;
    Tn = T(t);
    B = zeros(NDiv,NDiv);
    d = zeros(NDiv,1);     
    
    cinj1 = cinj(t);

    D = (D_ref/q_ref)*(0.5*(inj(t)+qp(t)));
    qf = (qf_ref/q_ref)*(0.5*(inj(t)+qp(t)));
    
    for m=1:NDiv
        Fx = (Vp*D)/(dx^2);
        cm = (Vp/dt);
        if m>1
            B(m,m-1) = (Fx+qf);
        end
        if m<NDiv
            B(m,m+1) = Fx;
        end
        if m==1
            B(m,m) = B(m,m)-Fx;
            d(m) = d(m) - (inj(t)*cinj1);
        end
        if m==NDiv
            B(m,m) = B(m,m) - Fx + qf;
            B(m,m) = B(m,m) + (-qp(t));
        end
        
        B(m,m) = B(m,m) - ((2*Fx)+qf+cm);
        d(m) = d(m) - (cm*C(m));
    end
        
    cwf_prod(t,1) = C(m);
    C = B\d;
end