PDEToolbox - Unsuitable initial guess U0 (default: U0=0)

Hello everyone.

I am trying to model a reaction of ammonia cracking ocurring in a Plug Flow (Cylindrical) Reactor with the geometry in the annex figure AQ4.jpg. For that end, I created the code below:

%1. Reactor Configuration [m]
L_aq4 = 1.5; %Reactor length
D_aq4 = 0.1; %Reactor diameter

%2. Miscellaneous Parameters
Diff_H2 = (5*10^-9)*800^1.682; %Average diffusion coefficient
lambda_gas = 6.09*10^-2; %Average thermal condutivity of the gas

%3. Creating the PDE Model
N = 5;
model = createpde(N);
% 1 = Total matter balance for the gases
% 2 = Matter balance for species NH3
% 3 = Matter balance for species N2
% 4 = Matter balance for species H2
% 5 = Heat balance

        %3.2.1. Creating the PDE Model - Geometry (Cylinder)
        gm = multicylinder (D_aq4/2,L_aq4);
        model.Geometry = gm;
        mesh_HminHmax = generateMesh(model,'Hmax',0.5,'Hmin',D_aq4);
        pdegplot(model,'FaceLabels','on','FaceAlpha',0.5);
        
        %3.2.2. Setting Boundary Conditions (Initial conditions not needed
        %for steady-state problems)
        C0_NH3 = 160; %[mol.m^-3]
        T0 = 300; %[K]
        Tw = 800; %[K]
        h1 = ones(N,1);
        h1d = diag(h1);
        h3 = [0 0 0 0 1];
        h3d = diag(h3);
        r1 = [C0_NH3 C0_NH3 0 0 T0];
        r3 = [0 0 0 0 Tw];
        applyBoundaryCondition(model,'dirichlet','Face',1,'h',h1d,'r',r1);
        applyBoundaryCondition(model,'dirichlet','Face',3,'h',h3d,'r',r3);
        
        %3.2.3. Setting PDE Coefficients (Elliptic equation type)
        c = [0 Diff_H2 Diff_H2 Diff_H2 lambda_gas]';
        specifyCoefficients(model,'m',0,'d',0,'c',c,'a',0,'f',@fcoeffunction);
        result = solvepde(model);
        pdeplot3D(model,'ColorMapData',result.NodalSolution)
        
        %4. Functions
function f = fcoeffunction(location,state)
nr = length(location.x);
f = zeros(5,nr);
R = 8.314; %Gas constant
ug = 1.5; %Gas velocity
cpavg_g = 3.727*10^1; %Average heat capacity of the gas
dHr_cr = 1.611*10^5; %Heat of cracking reaction
k_cr = 8.4*10^8; %Pre-exponential factor of cracking reaction
E_cr = 1.611*10^5; %Activation energy of cracking reaction
K_cr = k_cr.*exp(-E_cr./(R.*state.u(5,:))); %Arrhenius "constant" of cracking reaction
r_cr = K_cr.*state.u(2,:); %Reaction rate
f(1,:) = -ug.*state.uz(1,:) + 2.*r_cr;
f(2,:) = ug.*state.uz(2,:) + 2.*r_cr;
f(3,:) = ug.*state.uz(3,:) - r_cr;
f(4,:) = ug.*state.uz(4,:) - (3.*r_cr);
f(5,:) = (ug.*state.u(1,:).*cpavg_g).*state.uz(5,:) + dHr_cr.*r_cr;
end

The system of PDEs consist in the following equations:

1) ug*du(1)/dz = 2*r_cr

2) ug*du(2)/dz = -2*r_cr + Diff_H2*div(grad(u(2)))

3) ug*du(3)/dz = r_cr + Diff_H2*div(grad(u(3)))

4) ug*du(4)/dz = 3*r_cr + Diff_H2*div(grad(u(4)))

5) ug*cpavg_g*u(1)*du(5)/dz = -dHr_cr*r_cr + lambda*div(grad(u(5)))

For the boundary conditions:

a) For Face 1 (entrance of reactor): u(1) = u(2) = 160; u(3) = u(4) = 0; u(5) = 300;

b) For Face 3 (walls): u(5) = 800. Other variables are not bounded.

c) For Face 2 (exit): No boundary conditions.

But when I try to run the code, I get the following errors:

Warning: Matrix is singular to working precision.

> In pde.EquationModel/solveStationaryNonlinear (line 21)

In pde.PDEModel/solvepde (line 77)

In Codigo_AQ4 (line 55)

Error using pde.EquationModel/solveStationaryNonlinear (line 27)

Unsuitable initial guess U0 (default: U0=0).

By reading the documentation, the error "Unsuitable initial guess U0" means the initial guess is NaN or Inf. Is there any way to change the default U0 = 0 to prevent this? For what I understand, since I specified the boundary conditions at the entrace (Face 1), I would have specified U0, which, in this case, is not zero. Is there anything I am overlooking? Thanks in advance for the help.