# How to solve a system of six coupled PDE with step function

i'm trying to solve a coupled system of six partial differential equations with NDSolve. The system is as follows

but i don't have idea how to use the conditionals. Here $H_{\sigma}$ is a continuous mollifier of the step function satisfying

The parameter $\sigma > 0$ fixes the thickness of the transition between $H_\sigma(s) = 1$ and $H_\sigma(s) = 0$

My code is below

H[s_, \[Sigma]_] := Which[s <= 0, 1, s > \[Sigma], 0, True, 1 - s/\[Sigma]]

solucion[t_] =
With[{gaman = 0.746, gamat = 0.97, psi0 = 0.75, psi = 0.655,
psin = 0.6, psit = 0.8, deltat = 0.03, deltan = 0.1, mun = 0.1,
mut = 0.05, nu = 0.000016, pin = 6000000,
pit = 3000000, \[Tau] = 0.005, fo = 0.25, fg = 0.16, betan = 1.3,
betat = 1.2, Kan = 0.1, Kat = 0.3, \[Kappa] = 0.00000734, D = 1.0},

NDSolve[{
fin'[t] == oxi[t]*gaman*fin[t]*H[s, \[Sigma]] (psi - psin) - ((1 -
oxi[t]))*deltan*fin[t],

fit'[t] == glu[t]*gamat*fit[t]*H[s, \[Sigma]] (psi - psit) - ((1 -
glu[t]))*deltat*fit[t],

mat'[t] == mun*fin[t] + mun*fit[t] - nu*enz[t]*mat[t],

enz'[t] == pin*fin[t] + pit*fit[t],

oxi'[t] == -betan*fin[t]*oxi[t] + fo,

glu'[t] == -betat*fit[t]*glu[t] + fg,

fin[0] == 0.45, fit[0] == 0.005, mat[0] == 0.2, enz[0] == 0.2,
oxi[0] == 0.25, glu[0] == 0.16}, {fin[t], fit[t], mat[t], enz[t],
oxi[t], glu[t]}, {t, 0, 80}] /. {s -> 0, \[Sigma] -> 1}][[1, All, 2]]

Plot[solucion[t], {t, 0, 80}, PlotRange -> {{0, 80}, {0, .55}}]


but the above code does not correctly plot the expected results.

My explicit question is how to properly use conditionals inside NDSolve for $H_{\sigma}(s)$ for those three cases??

• Are you sure the definition of $H_{\sigma}(s)$ is correct? Its value is $1$ as long as $s\neq0$? Just try Plot3D[Evaluate@H[x, y], {x, -5, 5}, {y, -5, 5}] and you'll know what I mean. – xzczd Oct 25 '17 at 3:04
• What's the expected result? I think NDSolve has solved the given equation system correctly. – xzczd Oct 25 '17 at 5:21
• 1. The warning is simply because you don't place ReplaceAll (/.) in the correct place, but it doesn't cause mistake in the final result. If you don't want to see the warning, place the /. {s -> 0, \[Sigma] -> 1} inside NDSolve. 2. The symbol in your code is quite different from those in the picture, coding equations in this way is not a good idea because this makes error checking painful and inefficient. 3. The 4th equation is apparently different from the one in the picture. 4. You'd better add the Python code in your question if it's not too long. – xzczd Oct 26 '17 at 3:57