Errors obtained in the numerical solution of the 1D convection-diffusion equation

Question

I have been trying to solve the 1D convection-diffusion equation

In implementing the above, I have tried to follow closely the method described here: https://uk.mathworks.com/support/search.html/answers/391834-hi-i-m-trying-to-solve-the-heat-eq-using-the-explicit-and-implicit-methods-and-i-m-having-trouble-s.html?fq[]=asset_type_name:answer&fq[]=category:support/pde-solve5134&page=1

This is my code so far:

clc; clear; close all
 
%Input parameters
Ai     = 0; %Initial concentration in the column 
Ao    = 1;      %Influent concentration 
L     = 0.08;       %Length of column: m
nx    = 40;           %number of spatial gridpoints (i.e., no. of cells)
dx    = L/nx;           %Length step size [m]
T     = 20/24;      %End time [days]... that is 20hrs
nt    = 100;            %shifts
dt    = T/nt;       %Time step size [days] (Always check stability!)
Vel   = dx/dt;      %Velocity in each cell [m/day]
alpha = 0.002;      %Dispersivity [m]
De    = alpha*Vel;   % Dispersion coeff. [m2/day]
 
%Gridblocks
x    = 0:dx:L;
t    = 0:dt:T;
 
%Initial and boundary conditions
g1     = @(t) Ao;   
g2     = @(t) 0;
A       =  zeros(nx+1, nt+1);
A(:,1) =  Ai;     
A(1,:) =  g1(t);   

% Dimensionless parameters 
beta      = Vel*dt/dx;
gamma = De*dt/(dx^2);
 
% Implementation of the explicit method
for j= 1:nt-1      % Time Loop
    for i= 2:nx-1  % Space Loop
        A(i,j+1) = (A(i-1,j))*(beta + gamma)...
            + A(i,j)*(1-2*gamma-beta) + A(i+1,j)*(gamma);
    end
    
    % Insert boundary conditions for i = 1 and i = N
    A(2,j+1) =  A(1,j)*(beta + gamma) + A(2,j)*(1-2*gamma-beta) + A(3,j)*(gamma);
    A(nx,j+1) =  A(nx-1,j)*(beta + 2*gamma) + A(nx,j)*(1-2*gamma-Vel*beta)+ (2*gamma*dx*g2(t));
end
 
figure
plot(t, A(end,:), 'r*', 'MarkerSize', 2)
title('A Vs time profile (Using FDM)')
xlabel('t'),ylabel('A')

When I run the code, I obtain negative values of the dependent variable, however, comparison with MATLAB’s pdepe (see plot) gives a different profile, which I believe is the correct (matches with the analytical solution) result.

I am not sure what exactly I could be doing wrong here, although it seems that my major shortcoming is in the implementation of the initial and boundary conditions. I am thus hoping that someone may be able to help me out or point me in the right direction. Thank you in advance.