function dcdt = APBModel(t, c, grid, param)

% Initialize
    dcdt = zeros(grid.Ne, 1);
    ca = zeros(grid.N, 1);
    ga = c(grid.N-1:grid.Ne,1);

% Calculate BCs
    ca(1) = 1.0;
    ca(2:grid.N-1) = c(1:grid.N-2);

    isum = 0.0;
    for i=1:grid.N-1
        isum = isum + grid.A(grid.N,i)*ca(i);
    end
    ca(grid.N) = - isum/grid.A(grid.N,grid.N);


% Calculate the derivatives

    d2cdz2 = grid.B*ca;
    d1cdz1 = grid.A*ca;

% Model equations for interior points
for i=2:grid.N-1
    dcdt(i-1) = param.Pe*d2cdz2(i) - d1cdz1(i) - param.Film*(ca(i)-ga(i)/param.KcCa0/(1.0-ga(i)));
end


for i=1:grid.N
    dcdt(grid.N-2+i) = param.Gk*(ca(i)-ga(i)/param.KcCa0/(1.0-ga(i)));
end

end 

