Modeling and simulation of convection and diffusion is certainly possible to solve in Matlab with the FEA Toolbox, as shown in the model example below:
% Set up 1D domain from 0..1 with 20 elements.
fea.sdim = { 'x' };
fea.grid = linegrid( 20, 0, 1);
% Add covection and diffusion physics mode.
fea = addphys( fea, @convectiondiffusion, {'C'} );
% Define diffusion coefficient.
fea.phys.cd.eqn.coef{2,end} = {'d_coef'};
fea.expr = { 'c', {'1.23'} ;
'd_coef', {'(1+c)^2/((1+c)^2+1)'} };
% Use c = -1 on right boundary, and insulation
% flux boundary conditions on the left.
fea.phys.cd.bdr.sel = [ 1 3 ];
fea.phys.cd.bdr.coef{1,end}{1} = -1;
% Check, parse, and solve problem
% with initial condition 'C=2*x'.
fea = parsephys( fea );
fea = parseprob( fea );
[fea.sol.u,tlist] = ...
solvetime( fea, 'dt', 0.1, 'tmax', 1, 'init', {'2*x'} );
% Alternatively, solvestat can be used for stationary problems.
% Postprocessing.
isol = length(tlist);
postplot( fea, 'surfexpr', 'C', 'solnum', isol )
title( ['Solution at time, t = ',num2str(tlist(isol))] )
ylabel( 'Concentration, C' )
