# PDES+ODEs with ModelingToolkit

**URL:** <https://discourse.julialang.org/t/pdes-odes-with-modelingtoolkit/134481>\
**Category:** New to Julia\
**Created:** [December 10, 2025, 7:54pm UTC](https://discourse.julialang.org/t/pdes-odes-with-modelingtoolkit/134481 "2025-12-10T19:54:38Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![BPMD](https://avatars.discourse-cdn.com/v4/letter/b/f08c70/32.png) [@BPMD](https://discourse.julialang.org/u/BPMD)\
**Post date:** [December 10, 2025, 7:54pm UTC](https://discourse.julialang.org/t/pdes-odes-with-modelingtoolkit/134481/1 "2025-12-10T19:54:38Z")

</div>

I am simulating an electrochemical for Li+ purification using membranes. The problem includes several PDEs modeling the diffusion and migration in membrane coupled with ODEs representing the balances at chambers. The code currently is: using ModelingToolkit, DifferentialEquations, MethodOfLines, PDEBase, DomainSets  
using SciMLBase, LinearAlgebra, Plots, Sundials

# ----------------------

# 1. Physical and System Constants

# ----------------------

@parameters F D\_Li D\_OH RT z\_Li z\_OH k\_ext V\_anodo V\_catodo A  
@parameters F\_LiCl\_rec LiCl\_feed\_rec F\_LiOH\_rec Li\_OH\_rec F\_Wat\_Cat kf kr dPhi L  
@parameters C\_LiAIn C\_LiCIn C\_OHAIn C\_OHCIn

# ----------------------

# 1. Physical and System Constants (numeric)

# ----------------------

const F\_val = 96485.3329 # Faraday constant (C/mol)  
const R\_gas\_val = 8.3145 # Ideal gas constant (J/mol/K)  
const Temp\_val = 341.15 # Temperature (K)  
const D\_Li\_val = 9.0e-10 # Diffusion coefficient of Li+ (m^2/s)  
const D\_OH\_val = 2.7e-10 # Diffusion coefficient of OH- (m^2/s)  
const z\_Li\_val = 1 # Charge of Li+ ion  
const z\_OH\_val = -1 # Charge of OH- ion  
const L\_val = 1.5e-3 # Length of the domain (m)  
const N\_val = 20 # Number of spatial discretization points  
const k\_ext\_val = 2.0 # External mass transfer coefficient (m/s)

# ----------------------

# 2. Process Input/Flow Constants (numeric)

# ----------------------

const V\_anodo\_val = 15.0  
const V\_catodo\_val = 7.0  
const A\_val = 1.5  
const kf\_val = 1.0  
const kr\_val = 0.97  
const F\_LiCl\_rec\_val = 192/3600  
const F\_LiCl\_fre\_val = 0.15/3600  
const F\_Wat\_And\_val = 0.17/3600  
const F\_LiOH\_rec\_val = 192/3600  
const F\_Wat\_Cat\_val = 0.60/3600  
const LiCl\_feed\_rec\_val = 3.2  
const LiCl\_feed\_fre\_val = 11.83  
const Li\_OH\_rec\_val = 1.98

# ----------------------

# 3. Boundary/Initial Concentrations (numeric)

# ----------------------

const C\_LiAIn\_val = 0.25 # Li+ concentration at Anode Inlet  
const C\_LiCIn\_val = 0.15 # Li+ concentration at Cathode Inlet  
const C\_OHAIn\_val = 1e-4 # OH- concentration at Anode Inlet  
const C\_OHCIn\_val = 0.15 # OH- concentration at Cathode Inlet

# ----------------------

# 4. Derived Constants (numeric)

# ----------------------

const dx\_val = L\_val/(N\_val-1)  
const RT\_val = R\_gas\_val \* Temp\_val  
const RT\_over\_F\_val = RT\_val / F\_val

# Electric potentials

const Phi\_AIn\_val = RT\_over\_F\_val \* log(C\_LiAIn\_val / C\_LiCIn\_val)  
const Phi\_CIn\_val = 0.0

# Gradient of potential along x (numeric)

const dPhi\_dx\_val = (Phi\_CIn\_val - Phi\_AIn\_val)/L\_val

# Helper functions for Initial Conditions

C\_Li\_init\_func(x\_val) = C\_LiAIn + (C\_LiCIn-C\_LiAIn)\*x\_val/L  
C\_OH\_init\_func(x\_val) = C\_OHAIn + (C\_OHCIn-C\_OHAIn)\*x\_val/L

# ----------------------

# 5. Symbolic PDEs Definition

# ----------------------

@variables x  
@independent\_variables t

# PDE variables

@variables C\_Li(..) C\_OH(..)

# ODE variables

@variables C\_LiA(..) C\_LiC(..) C\_LiOHC(..) C\_OHC(..)

# Fluxes and reaction rate (functions of time)

@variables flux\_Li\_anod(..) flux\_Li\_cathod(..) flux\_OH\_cathod(..) reaction\_rate\_cathode(..)

Dt = Differential(t)  
Dx = Differential(x)  
Dxx = Differential(x)^2

# PDEs

pde\_eqs = [  
Dt(C\_Li(x,t)) ~ D\_Li_Dx(D\_x(C\_Li(x,t))) + D\_Li_(F/RT)_dPhi_Dx(C\_Li(x,t)),  
Dt(C\_OH(x,t)) ~ D\_OH_Dx(D\_x(C\_OH(x,t))) - D\_OH_(F/RT)_dPhi_Dx(C\_OH(x,t))  
]

# Flux/reaction equations

anode\_flux\_eq\_Li = flux\_Li\_anod(t) ~ -D\_Li \* Dx(C\_Li(0, t)) - (z\_Li \* D\_Li \* F / RT) \* C\_Li(0, t) \* dPhi  
cathod\_flux\_eq\_Li = flux\_Li\_cathod(t) ~ -D\_Li \* Dx(C\_Li(L, t)) - (z\_Li \* D\_Li \* F / RT) \* C\_Li(L, t) \* dPhi  
cathod\_flux\_eq\_OH = flux\_OH\_cathod(t) ~ -D\_OH \* Dx(C\_OH(L, t)) - (z\_OH \* D\_OH \* F / RT) \* C\_OH(L, t) \* dPhi  
reaction\_rate\_eq = reaction\_rate\_cathode(t) ~ kf \* C\_LiC(t) \* C\_OHC(t) - kr \* C\_LiOHC(t)

# ODEs

ode\_eqs = [  
Dt(C\_LiA(t)) ~ (F\_LiCl\_rec/V\_anodo)_(LiCl\_feed\_rec - C\_LiA(t)) - (A/V\_anodo) \* flux\_Li\_anod(t),  
Dt(C\_LiC(t)) ~ (F\_LiOH\_rec/V\_catodo)_(Li\_OH\_rec - C\_LiC(t)) + (A/V\_catodo) \* flux\_Li\_cathod(t) - reaction\_rate\_cathode(t),  
Dt(C\_LiOHC(t)) ~ (F\_LiOH\_rec/V\_catodo) \* Li\_OH\_rec - ((F\_LiOH\_rec + F\_Wat\_Cat) / V\_catodo) \* C\_LiOHC(t) + reaction\_rate\_cathode(t),  
Dt(C\_OHC(t)) ~ (F\_LiOH\_rec/V\_catodo) \* Li\_OH\_rec - ((F\_LiOH\_rec + F\_Wat\_Cat) / V\_catodo) \* C\_OHC(t) + (A/V\_catodo) \* flux\_OH\_cathod(t) - reaction\_rate\_cathode(t)  
]

# Combine all equations

eqs = vcat(pde\_eqs, [anode\_flux\_eq\_Li, cathod\_flux\_eq\_Li, cathod\_flux\_eq\_OH, reaction\_rate\_eq], ode\_eqs)

bcs = [  
# Initial Conditions (t=0) for PDEs  
C\_Li(x, 0) ~ C\_Li\_init\_func(x),  
C\_OH(x, 0) ~ C\_OH\_init\_func(x),

```
# Initial Conditions for the ODE variables
C_LiA(0) ~ C_LiAIn,
C_LiC(0) ~ C_LiCIn,
C_OHC(0) ~ C_OHCIn,
C_LiOHC(0) ~ 0.0,

# Boundary Flux Conditions (x=0, Anode side)
Dx(C_Li(0,t)) + (z_Li*F/RT)*C_Li(0,t)*dPhi ~ -k_ext/D_Li*(C_Li(0,t) - C_LiA(t)),
Dx(C_OH(0,t)) + (z_OH*F/RT)*C_OH(0,t)*dPhi ~ -k_ext/D_OH*(C_OH(0,t) - C_OHAIn),

# Boundary Flux Conditions (x=L, Cathode side)
Dx(C_Li(L,t)) + (z_Li*F/RT)*C_Li(L,t)*dPhi ~ -k_ext/D_Li*(C_Li(L,t) - C_LiC(t)),
Dx(C_OH(L,t)) + (z_OH*F/RT)*C_OH(L,t)*dPhi ~ -k_ext/D_OH*(C_OH(L,t) - C_OHC(t))

```

]

# ----------------------

# 6. System Construction

# ----------------------

domains = [x ∈ Interval(0.0, L\_val), t ∈ Interval(0.0, 100.0)]

# Collect all parameters needed in the defaults dictionary

all\_params = Dict(  
F =\> F\_val,  
D\_Li =\> D\_Li\_val,  
D\_OH =\> D\_OH\_val,  
RT =\> RT\_val,  
z\_Li =\> z\_Li\_val,  
z\_OH =\> z\_OH\_val,  
k\_ext =\> k\_ext\_val,  
V\_anodo =\> V\_anodo\_val,  
V\_catodo =\> V\_catodo\_val,  
A =\> A\_val,  
F\_LiCl\_rec =\> F\_LiCl\_rec\_val,  
LiCl\_feed\_rec =\> LiCl\_feed\_rec\_val,  
F\_LiOH\_rec =\> F\_LiOH\_rec\_val,  
Li\_OH\_rec =\> Li\_OH\_rec\_val,  
F\_Wat\_Cat =\> F\_Wat\_Cat\_val,  
C\_LiAIn =\> C\_LiAIn\_val,  
C\_LiCIn =\> C\_LiCIn\_val,  
C\_OHAIn =\> C\_OHAIn\_val,  
C\_OHCIn =\> C\_OHCIn\_val,  
kf =\> kf\_val,  
kr =\> kr\_val,  
dPhi =\> dPhi\_dx\_val,  
L =\> L\_val,  
)

@named pde\_system = PDESystem(  
eqs,  
bcs,  
domains,  
[x, t],  
[C\_Li(x,t), C\_OH(x,t), C\_LiA(t), C\_LiC(t), C\_OHC(t), C\_LiOHC(t),  
flux\_Li\_anod(t), flux\_Li\_cathod(t), flux\_OH\_cathod(t), reaction\_rate\_cathode(t)],  
collect(keys(all\_params));  
defaults=all\_params  
)

# ----------------------

# 7. Discretization and Solution

# ----------------------

@show dx\_val  
xspan = 0.0:Float64(dx\_val):Float64(L\_val)  
discretization = MOLFiniteDifference([x =\> dx\_val], t; advection\_scheme=UpwindScheme())  
prob\_ode = discretize(pde\_system, discretization)

# Solve the problem

sol = solve(prob\_ode, TRBDF2(), saveat=1.0, reltol=1e-6, abstol=1e-8)

# ----------------------

# 8. Results Visualization

# ----------------------

# Extract discrete x values

x\_vals = range(0, stop=L, length=N)

# Extract concentration matrices

C\_Li\_mat = sol[C\_Li(x, t)]  
C\_OH\_mat = sol[C\_OH(x, t)]

# Target times for plotting

target\_times = [0.0, 25.0, 50.0, 75.0, 100.0]

# Create plots for each time

p = plot(layout=(length(target\_times), 1), size=(800, 1000))

for (i, t\_target) in enumerate(target\_times)  
# Find closest time index  
t\_idx = argmin(abs.(sol.t .- t\_target))

```
# Extract profiles at that time
Li_profile = C_Li_mat[:, t_idx]
OH_profile = C_OH_mat[:, t_idx]

# Plot both species
plot!(
    p[i],
    x_vals, Li_profile,
    label="Li⁺",
    lw=2,
    color=:blue,
    title="t = $(t_target) s",
    legend=:topright
)
plot!(
    p[i],
    x_vals, OH_profile,
    label="OH⁻",
    lw=2,
    color=:red,
    legend=:topright
)

# Format axes
xlabel!(p[i], "Position (m)")
ylabel!(p[i], "Concentration (mol/m³)")

```

end

# Save the figure

savefig(p, “concentration\_profiles.png”)  
but I got an error: ERROR: LoadError: ArgumentError: invalid index: nothing of type Nothing Stacktrace: [1] to\_index(i::Nothing) @ Base .\indices.jl:315 [2] to\_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing) @ Base .\indices.jl:292 [3] to\_indices @ .\indices.jl:368 [inlined] [4] to\_indices @ .\indices.jl:360 [inlined] [5] getindex @ .\abstractarray.jl:1312 [inlined] [6] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real}) The error is in prob\_ode = discretize(pde\_system, discretization) line. Previously I solved the problem including only the membrane. But when the balances are included the BC’s of the PDEs are intelinked with them. If someone has ideas how to handle the problem, I will be very grateful.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 10, 2025, 8:14pm UTC](https://discourse.julialang.org/t/pdes-odes-with-modelingtoolkit/134481/2 "2025-12-10T20:14:55Z")

</div>

**Small remark:**  
If you share code, please enclose it in triple backticks, like ```

That makes it much easier to read your code, but also to copy-and-paste it for people who want to reproduce your problem.

---

<div class="post-metadata">

**Author:** ![BPMD](https://avatars.discourse-cdn.com/v4/letter/b/f08c70/32.png) [@BPMD](https://discourse.julialang.org/u/BPMD)\
**Post date:** [December 10, 2025, 10:14pm UTC](https://discourse.julialang.org/t/pdes-odes-with-modelingtoolkit/134481/3 "2025-12-10T22:14:25Z")

</div>

```using
using SciMLBase, LinearAlgebra, Plots, Sundials

# ----------------------
# 1. Physical and System Constants
# ----------------------
@parameters F D_Li D_OH RT z_Li z_OH k_ext V_anodo V_catodo A 
@parameters F_LiCl_rec LiCl_feed_rec F_LiOH_rec Li_OH_rec F_Wat_Cat kf kr dPhi L
@parameters C_LiAIn C_LiCIn C_OHAIn C_OHCIn

# ----------------------
# 1. Physical and System Constants (numeric)
# ----------------------
const F_val = 96485.3329 # Faraday constant (C/mol)
const R_gas_val = 8.3145 # Ideal gas constant (J/mol/K)
const Temp_val = 341.15 # Temperature (K)
const D_Li_val = 9.0e-10 # Diffusion coefficient of Li+ (m^2/s)
const D_OH_val = 2.7e-10 # Diffusion coefficient of OH- (m^2/s)
const z_Li_val = 1 # Charge of Li+ ion
const z_OH_val = -1 # Charge of OH- ion
const L_val = 1.5e-3 # Length of the domain (m)
const N_val = 20 # Number of spatial discretization points
const k_ext_val = 2.0 # External mass transfer coefficient (m/s)

# ----------------------
# 2. Process Input/Flow Constants (numeric)
# ----------------------
const V_anodo_val = 15.0
const V_catodo_val = 7.0
const A_val = 1.5
const kf_val = 1.0
const kr_val = 0.97
const F_LiCl_rec_val = 192/3600
const F_LiCl_fre_val = 0.15/3600
const F_Wat_And_val = 0.17/3600
const F_LiOH_rec_val = 192/3600
const F_Wat_Cat_val = 0.60/3600
const LiCl_feed_rec_val = 3.2
const LiCl_feed_fre_val = 11.83
const Li_OH_rec_val = 1.98

# ----------------------
# 3. Boundary/Initial Concentrations (numeric)
# ----------------------
const C_LiAIn_val = 0.25 # Li+ concentration at Anode Inlet
const C_LiCIn_val = 0.15 # Li+ concentration at Cathode Inlet
const C_OHAIn_val = 1e-4 # OH- concentration at Anode Inlet
const C_OHCIn_val = 0.15 # OH- concentration at Cathode Inlet

# ----------------------
# 4. Derived Constants (numeric)
# ----------------------
const dx_val = L_val/(N_val-1)
const RT_val = R_gas_val * Temp_val
const RT_over_F_val = RT_val / F_val

# Electric potentials
const Phi_AIn_val = RT_over_F_val * log(C_LiAIn_val / C_LiCIn_val)
const Phi_CIn_val = 0.0

# Gradient of potential along x (numeric)
const dPhi_dx_val = (Phi_CIn_val - Phi_AIn_val)/L_val

# Helper functions for Initial Conditions
C_Li_init_func(x_val) = C_LiAIn + (C_LiCIn-C_LiAIn)*x_val/L
C_OH_init_func(x_val) = C_OHAIn + (C_OHCIn-C_OHAIn)*x_val/L

# ----------------------
# 5. Symbolic PDEs Definition
# ----------------------

@variables x
@independent_variables t
# PDE variables
@variables C_Li(..) C_OH(..)

# ODE variables
@variables C_LiA(..) C_LiC(..) C_LiOHC(..) C_OHC(..)

# Fluxes and reaction rate (functions of time)
@variables flux_Li_anod(..) flux_Li_cathod(..) flux_OH_cathod(..) reaction_rate_cathode(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# PDEs
pde_eqs = [
    Dt(C_Li(x,t)) ~ D_Li*Dx(D_x(C_Li(x,t))) + D_Li*(F/RT)*dPhi*Dx(C_Li(x,t)),
    Dt(C_OH(x,t)) ~ D_OH*Dx(D_x(C_OH(x,t))) - D_OH*(F/RT)*dPhi*Dx(C_OH(x,t))
]

# Flux/reaction equations
anode_flux_eq_Li = flux_Li_anod(t) ~ -D_Li * Dx(C_Li(0, t)) - (z_Li * D_Li * F / RT) * C_Li(0, t) * dPhi
cathod_flux_eq_Li = flux_Li_cathod(t) ~ -D_Li * Dx(C_Li(L, t)) - (z_Li * D_Li * F / RT) * C_Li(L, t) * dPhi
cathod_flux_eq_OH = flux_OH_cathod(t) ~ -D_OH * Dx(C_OH(L, t)) - (z_OH * D_OH * F / RT) * C_OH(L, t) * dPhi
reaction_rate_eq = reaction_rate_cathode(t) ~ kf * C_LiC(t) * C_OHC(t) - kr * C_LiOHC(t)

# ODEs
ode_eqs = [
    Dt(C_LiA(t)) ~ (F_LiCl_rec/V_anodo)*(LiCl_feed_rec - C_LiA(t)) - (A/V_anodo) * flux_Li_anod(t),
    Dt(C_LiC(t)) ~ (F_LiOH_rec/V_catodo)*(Li_OH_rec - C_LiC(t)) + (A/V_catodo) * flux_Li_cathod(t) - reaction_rate_cathode(t),
    Dt(C_LiOHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_LiOHC(t) + reaction_rate_cathode(t),
    Dt(C_OHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_OHC(t) + (A/V_catodo) * flux_OH_cathod(t) - reaction_rate_cathode(t)
]

# Combine all equations
eqs = vcat(pde_eqs, [anode_flux_eq_Li, cathod_flux_eq_Li, cathod_flux_eq_OH, reaction_rate_eq], ode_eqs)

bcs = [
    # Initial Conditions (t=0) for PDEs
    C_Li(x, 0) ~ C_Li_init_func(x),
    C_OH(x, 0) ~ C_OH_init_func(x),
    
    # Initial Conditions for the ODE variables
    C_LiA(0) ~ C_LiAIn,
    C_LiC(0) ~ C_LiCIn,
    C_OHC(0) ~ C_OHCIn,
    C_LiOHC(0) ~ 0.0,
    
    # Boundary Flux Conditions (x=0, Anode side)
    Dx(C_Li(0,t)) + (z_Li*F/RT)*C_Li(0,t)*dPhi ~ -k_ext/D_Li*(C_Li(0,t) - C_LiA(t)),
    Dx(C_OH(0,t)) + (z_OH*F/RT)*C_OH(0,t)*dPhi ~ -k_ext/D_OH*(C_OH(0,t) - C_OHAIn),
    
    # Boundary Flux Conditions (x=L, Cathode side)
    Dx(C_Li(L,t)) + (z_Li*F/RT)*C_Li(L,t)*dPhi ~ -k_ext/D_Li*(C_Li(L,t) - C_LiC(t)),
    Dx(C_OH(L,t)) + (z_OH*F/RT)*C_OH(L,t)*dPhi ~ -k_ext/D_OH*(C_OH(L,t) - C_OHC(t))
]
# ----------------------
# 6. System Construction
# ----------------------
domains = [x ∈ Interval(0.0, L_val), t ∈ Interval(0.0, 100.0)]

# Collect all parameters needed in the defaults dictionary
all_params = Dict(
    F => F_val,
    D_Li => D_Li_val,
    D_OH => D_OH_val,
    RT => RT_val,
    z_Li => z_Li_val,
    z_OH => z_OH_val,
    k_ext => k_ext_val,
    V_anodo => V_anodo_val,
    V_catodo => V_catodo_val,
    A => A_val,
    F_LiCl_rec => F_LiCl_rec_val,
    LiCl_feed_rec => LiCl_feed_rec_val,
    F_LiOH_rec => F_LiOH_rec_val,
    Li_OH_rec => Li_OH_rec_val,
    F_Wat_Cat => F_Wat_Cat_val,
    C_LiAIn => C_LiAIn_val, 
    C_LiCIn => C_LiCIn_val,
    C_OHAIn => C_OHAIn_val, 
    C_OHCIn => C_OHCIn_val,
    kf => kf_val,
    kr => kr_val,
    dPhi => dPhi_dx_val,
    L => L_val,
)

@named pde_system = PDESystem(
    eqs,
    bcs,
    domains,
    [x, t],
    [C_Li(x,t), C_OH(x,t), C_LiA(t), C_LiC(t), C_OHC(t), C_LiOHC(t), 
    flux_Li_anod(t), flux_Li_cathod(t), flux_OH_cathod(t), reaction_rate_cathode(t)],
    collect(keys(all_params));
    defaults=all_params
)

# ----------------------
# 7. Discretization and Solution
# ----------------------
@show dx_val
xspan = 0.0:Float64(dx_val):Float64(L_val)
discretization = MOLFiniteDifference([x => dx_val], t; advection_scheme=UpwindScheme())
prob_ode = discretize(pde_system, discretization)

# Solve the problem
sol = solve(prob_ode, TRBDF2(), saveat=1.0, reltol=1e-6, abstol=1e-8)

# ----------------------
# 8. Results Visualization
# ----------------------
# Extract discrete x values
x_vals = range(0, stop=L, length=N)

# Extract concentration matrices
C_Li_mat = sol[C_Li(x, t)]
C_OH_mat = sol[C_OH(x, t)]

# Target times for plotting
target_times = [0.0, 25.0, 50.0, 75.0, 100.0]

# Create plots for each time
p = plot(layout=(length(target_times), 1), size=(800, 1000))

for (i, t_target) in enumerate(target_times)
    # Find closest time index
    t_idx = argmin(abs.(sol.t .- t_target))
    
    # Extract profiles at that time
    Li_profile = C_Li_mat[:, t_idx]
    OH_profile = C_OH_mat[:, t_idx]
    
    # Plot both species
    plot!(
        p[i],
        x_vals, Li_profile,
        label="Li⁺",
        lw=2,
        color=:blue,
        title="t = $(t_target) s",
        legend=:topright
    )
    plot!(
        p[i],
        x_vals, OH_profile,
        label="OH⁻",
        lw=2,
        color=:red,
        legend=:topright
    )
    
    # Format axes
    xlabel!(p[i], "Position (m)")
    ylabel!(p[i], "Concentration (mol/m³)")
end

# Save the figure
savefig(p, "concentration_profiles.png")

```