How to add multiple Y axis to the same graph

Hello, I have a working PDE that plots the concentration profile of CO2 released by roots in a soil column by depth, I am wanting to add a second axis so I can plot what the exponential decay function (beta) looks like on top of my concentration profile. What would be the most straightforward way to do this? I am very new to matlab so this is a bit of a learning curve for me. I have read a little about yyaxis functions but I am not entirely sure how to impliment this.

function diffusion_roots_pdepe

x = linspace(0, 1, 500); %divides soil depth of 1m into 500 points
t = linspace(0, 86400, 500); %defines time points where concentration is computed
m = 0; %1D slab (cartesian coordinates)



sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);
C = sol(:, :,1); %Solve for 1 variable across all time and space points

figure
plot(x, C(end,:), 'LineWidth', 2)
xlabel('Distance (m)')
ylabel('Concentration (mol m⁻³)')
title('Final Concentration Profile')
grid on

figure
h_line = plot(x, C(1,:), 'b-', 'LineWidth', 2);
xlabel('Depth (m)')
ylabel('Concentration (mol m^{-3})')
title('Soil gas concentration profile')
grid on

ylim([min(C(:)), max(C(:))])
xlim([0, 1])
h_title = title('');  

for i = 1:length(t)
    set(h_line, 'YData', C(i,:))
    h_title.String = sprintf('CO2 concentration — t = %.0f s', t(i));
    drawnow
end

function [c, f, s ] = pdefun(x, t, u, dudx)

    %c = storage term
    %f = flux term
    %s = source/sink   

%EFFECTIVE AIR DIFFUSIVITY

%BUCKINGHAM DIFFUSIVITY - SANDY SOILS

D_air = 1.6e-5; %Diffusion of CO2 in free air, m² s⁻¹
theta = 0.6; %Air filled porosity of soil 
Dp_Do = theta^2; %Buckingham relative soil diffusivity, where theta is gas filled porosity of soil
D_eff = Dp_Do * D_air; %Effective soil diffusivity scaled by Dp_Do

%MOLDRUP DIFFUSIVITY - HEAVY CLAY SOILS

epsilon = 0.12 %Air filled porosity
phi  = 0.53 %total soil porosity 
b = 11.4 %campbell soil water retention parameter

D_eff_clay_1 = D_air * (phi^2)  * (epsilon / phi)^(2 + 3/b);

%MILLINTON-QUIRK DIFFUSIVITY

D_eff_clay_2 = D_air * (epsilon^2) / (phi ^ (2/3));

%ROOT DISTRIBUTION

beta = 5.961; %root decay factor defining distribution with depth
R = exp(-beta * x); %Gale Grigal formula for modelling root density with depth

%ROOT EMISSION 


S0 = 1e-8; % mol m⁻³ s⁻¹ — root emission rate
S = S0*R; %Scales production by root system distribution

%MICROBIAL SINK

K_micro = 0; %Microbial consumption rate, mol m⁻³ s⁻¹ (Vmax)
Km = 1e-6; %Saturation point, when concentration is half Vmax, mol m⁻³
%This results in the rate of increase in consumption to slow down until
%it approaches saturation and slows down 

sink = K_micro * R * u / (Km + u); %Michaelis-Menten Kinetics term 
%for non linear consumption limited by saturation, leads to saturation of
%uptake by microbes at high concentration above Km


%PDE TERMS

c = 1; %Storage coefficient 
f = D_eff * dudx; %Fickian diffusion term scaled by effective soil diffusivity
s = S - sink; %Production - Consumption 

end

%STARTING CONDITIONS
    function U0 = icfun(x) 
        U0 = 0.016; %Initial concentration, atmospheric CO2
    end

%BOUNDARY CONDITIONS 

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
    %pl, ql, pr, qr left and right boundary coefficients respectively

    h = 0.00005; %Boundary exchange coefficient, m s⁻¹
    C_atm = 0.016; %Atmospheric concentration, mol m⁻³

    % Left boundary Robin condition
    pl = -h*(ul - C_atm); %Robin coefficient, where diffusion to atmosphere
    %is limited by h, a factor at the boundary restricting diffusion 
    ql = 1;  % Flux across the boundary
    % Right boundary Neumann condition
    pr = 0; %there is no flux across the boundary, so no flux term 
    qr = 1;%Flux term active at the boundary so as to not fix concentration
    %There is no flux across a Neumann boundary, this can lead to build up
    %of substrate

  
end

end

%structure of code was based on matlab documentation on how to use PDEPE 

Answer 1

[dpb]

Try the following modifications -- I commented out the second figure because in the end it duplicates the first. For display purposes, you might put a pause in that loop so you could watch the evolution in time or look up how to make a video by stringing frames together.

I made an assumption the exponential wanted is the variable R down in the pde function; for expediency here I simply copied over the code and turned it into an anonymous function; in your code in a better editing environment, you can make a separate public function or use one of the ways to pass additional variables to/from a pde solver. It probably needs a scale factor applied as well to turn it into an actual density; I didn't dig into the code enough to try to figure out where that may be.

yyaxis is really quite simple as shown -- just call it first to set which axis is the target; as you can see, there are then two YAxis objects in the axes object; the first is left, the second right. Used here to fix the ugly variable number of decimal places on the tick labels the default '%g' format does. You'll need to fix up the units in the label.

diffusion_roots_pdepe

function diffusion_roots_pdepe

% make a temporary anonymous function for plotting purposes for here...
beta = 5.961;           %root decay factor defining distribution with depth
fnR=@(b,x)exp(-b*x);    %Gale Grigal formula for modelling root density with depth

x = linspace(0, 1, 500); %divides soil depth of 1m into 500 points
t = linspace(0, 86400, 500); %defines time points where concentration is computed
m = 0; %1D slab (cartesian coordinates)

C=pdepe(m, @pdefun, @icfun, @bcfun, x, t);
% solution is C, no need for the temporary variable, it is just duplicated
%C = sol(:, :,1); %Solve for 1 variable across all time and space points

figure
yyaxis left                             % plot normal for concentration first
plot(x, C(end,:), 'LineWidth', 2)
xlabel('Distance (m)')
ylabel('Concentration (mol m⁻³)')
title('Final Concentration Profile')
hAx=gca;
hAx.YAxis(1).TickLabelFormat='%0.6f';
grid on
% now do the root density stuff...
R=fnR(beta,x);
yyaxis right
plot(x,R,'LineWidth',2)
ylabel('Root density (units(?) m⁻³)')
hAx.YAxis(2).TickLabelFormat='%0.2f';
%plot()
% ends up as duplicate of first figure, don't need for now...
% figure
% h_line = plot(x, C(1,:), 'b-', 'LineWidth', 2);
% xlabel('Depth (m)')
% ylabel('Concentration (mol m^{-3})')
% title('Soil gas concentration profile')
% grid on
% 
% ylim([min(C(:)), max(C(:))])
% xlim([0, 1])
% h_title = title('');  
% 
% for i = 1:length(t)
%     set(h_line, 'YData', C(i,:))
%     h_title.String = sprintf('CO2 concentration — t = %.0f s', t(i));
%     drawnow
% end

function [c, f, s ] = pdefun(x, t, u, dudx)

    %c = storage term
    %f = flux term
    %s = source/sink   

%EFFECTIVE AIR DIFFUSIVITY

%BUCKINGHAM DIFFUSIVITY - SANDY SOILS

D_air = 1.6e-5; %Diffusion of CO2 in free air, m² s⁻¹
theta = 0.6; %Air filled porosity of soil 
Dp_Do = theta^2; %Buckingham relative soil diffusivity, where theta is gas filled porosity of soil
D_eff = Dp_Do * D_air; %Effective soil diffusivity scaled by Dp_Do

%MOLDRUP DIFFUSIVITY - HEAVY CLAY SOILS

epsilon = 0.12; %Air filled porosity
phi  = 0.53;    %total soil porosity 
b = 11.4;       %campbell soil water retention parameter

D_eff_clay_1 = D_air * (phi^2)  * (epsilon / phi)^(2 + 3/b);

%MILLINTON-QUIRK DIFFUSIVITY

D_eff_clay_2 = D_air * (epsilon^2) / (phi ^ (2/3));

%ROOT DISTRIBUTION

beta = 5.961; %root decay factor defining distribution with depth
R = exp(-beta * x); %Gale Grigal formula for modelling root density with depth

%ROOT EMISSION 


S0 = 1e-8; % mol m⁻³ s⁻¹ — root emission rate
S = S0*R; %Scales production by root system distribution

%MICROBIAL SINK

K_micro = 0; %Microbial consumption rate, mol m⁻³ s⁻¹ (Vmax)
Km = 1e-6; %Saturation point, when concentration is half Vmax, mol m⁻³
%This results in the rate of increase in consumption to slow down until
%it approaches saturation and slows down 

sink = K_micro * R * u / (Km + u); %Michaelis-Menten Kinetics term 
%for non linear consumption limited by saturation, leads to saturation of
%uptake by microbes at high concentration above Km


%PDE TERMS

c = 1; %Storage coefficient 
f = D_eff * dudx; %Fickian diffusion term scaled by effective soil diffusivity
s = S - sink; %Production - Consumption 

end

%STARTING CONDITIONS
    function U0 = icfun(x) 
        U0 = 0.016; %Initial concentration, atmospheric CO2
    end

%BOUNDARY CONDITIONS 

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
    %pl, ql, pr, qr left and right boundary coefficients respectively

    h = 0.00005; %Boundary exchange coefficient, m s⁻¹
    C_atm = 0.016; %Atmospheric concentration, mol m⁻³

    % Left boundary Robin condition
    pl = -h*(ul - C_atm); %Robin coefficient, where diffusion to atmosphere
    %is limited by h, a factor at the boundary restricting diffusion 
    ql = 1;  % Flux across the boundary
    % Right boundary Neumann condition
    pr = 0; %there is no flux across the boundary, so no flux term 
    qr = 1;%Flux term active at the boundary so as to not fix concentration
    %There is no flux across a Neumann boundary, this can lead to build up
    %of substrate

  
end

end

%structure of code was based on matlab documentation on how to use PDEPE