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);
for t = 1:(length(T))
    dt = 1;
    B = zeros(NDiv,NDiv);
    d = zeros(NDiv,1);
    
    if t<=1000
        cinj1 = 250;
    elseif t>1000 && t<=2000
        cinj1 = 200;
    elseif t>2000 && t<=2500
        cinj1 = 350;
    elseif t>2500 && t<=3500
        cinj1 = 170;
    elseif t>3500 && t<=5000
        cinj1 = 200;
    elseif t>5000 && t<=7500
        cinj1 = 250;
    elseif t>7500 && t<=10000
        cinj1 = 350;
    elseif t>10000 && t<=15000
        cinj1 = 220;
    elseif t>15000 && t<=18000
        cinj1 = 300;
    elseif t>18000
        cinj1 = 265;
    end
%     cinj1 = cinj;

    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