How to plot parameters within bvp4c

Hello,

I'm solving for a reaction diffusion equation. I followed examples on how to plot the solution of differential equation and was able to do that. However, I'm not too sure how to plot the N_r, N_R, N_r2, N_rR, N_pR, N_R2, R_L relative to C_L or x.

function time_dependent_pdepe
clear all; close all; clc;

D_ij = 30; %Diffusion coefficient D (3.0*10^-7 cm^2/s -> 30 um^2/s)
L0 = 10; %c0 [nM]
x_f =40; %Length of domain [um]
maxt = 10; %Max simulation time [s]

m = 0; %Parameter corresponding to the symmetry of the problem
x = linspace(0,x_f,100); %xmesh
t = linspace(0,maxt,100); %tspan

sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
u = sol;

% Plotting
hold all
for n = linspace(1,length(t),10)
plot(x,sol(n,:),'LineWidth',2)

end

xlabel('Distance (\mum)')
ylabel('Concentration (nM)')
axis([0 x_f 0 L0])

function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)
D = D_ij;

%Rate constants
k_1 = 0.00193;
k_2 = 0.00255;
k_3 = 4.09;
d_1 = 0.007;
d_2 = 3.95*10^-5;
d_3 = 2.26;

%Equilibrium constants
K_1 = 3.63;
K_2 = 0.0155;
K_3 = 0.553;
K_4 = 9.01;
K_i = 0.139;

%Total Recptors
N_T = 1.7;

N_r = -(K_3.*K_4.*K_i.*(K_2 + u - ((8.*K_2.*u.^2.*N_T + 8.*K_2.*K_3.*u.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*u.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*u + 8.*K_2.*K_3.*K_i.*u.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*u.^2 + 4.*K_3.*u + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*u);
N_R = (u./K_1).*N_r;
N_r2 = (1./K_4).*(N_r).^2;
N_rR = (u./(2.*K_2.*K_4)).*(N_r).^2;
N_pR = (u./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
N_R2 = ((u.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*u.*N_r))+...
    ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*u.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
    (d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*u.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

% PDE
c = 1;
f = D_ij.*dudx;
s = R_L;
 end

function u0 = DiffusionICfun(x)
u0 = 0;
end

function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)

c0 = L0;

% At x = 0, c0 = L0
pl = ul-c0; 
ql = 0; 

% At x = L, dc/dx = 0
pr = 0;
qr = 1;

end

end

When I try plotting it, I get a blank plot. How would you go about plotting those values?

Thank you!