As per your previous question https://www.mathworks.com/matlabcentral/answers/368440-how-to-solve-fick-s-2nd-law-of-diffusion-equation
% Set up 1D domain from 0..3 with 20 elements.
fea.sdim = { 'x' };
fea.grid = linegrid( 20, 0, 3 );
% Add covection and diffusion physics mode.
fea = addphys( fea, @convectiondiffusion, {'u'} );
% Define diffusion coefficient.
fea.phys.cd.eqn.coef{2,end} = {'d_coef'};
fea.expr = { 'd_coef', {'0.00296*(1+49*u)^2/((1+49*u)^2+0.00164)'} };
% Use u = 0 on left boundary (x=0), and insulation
% flux boundary du/dx=0 conditions on the right (x=3).
fea.phys.cd.bdr.sel = [ 1 3 ];
fea.phys.cd.bdr.coef{1,end}{1} = 0;
% Check, parse, and solve problem
% with initial condition 'u=0.01'.
fea = parsephys( fea );
fea = parseprob( fea );
[fea.sol.u,tlist] = ...
solvetime( fea, 'tstep', 10, 'tmax', 100, 'init', {'0.01'} );
% Alternatively, solvestat can be used for stationary problems.
% Postprocessing.
isol = length(tlist);
postplot( fea, 'surfexpr', 'u', 'solnum', 1 )
hold on
postplot( fea, 'surfexpr', 'u', 'solnum', floor(isol/2) )
postplot( fea, 'surfexpr', 'u', 'solnum', isol )
axis normal
grid on
title( ['Solution at times, t = ',num2str(tlist([1 floor(isol/2) isol]))] )
ylabel( 'Concentration, u' )
