Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm solving some coupled PDEs (Eb1 and Ef1) and what I plot for Eb1 appears to be correct. However, for some reason, when I go to plot Ef1, I get nothing. MWE is below (beginning is constants and functions needed to solve PDEs, relevant part near end)

(*Constants needed*)
a = 8.2314*10^-7; omega = 3.0318*10^7; Do2 = 2*10^-9; po = 106;
ro = 235*10^-6; micron = 1*10^-6; qm = 10^-4;
mic = 1*10^-6; De = 5.5*10^-11; eo = 100; j = 0.12; kmn = 4.7;
kme = 2.1; con = (a*omega)/(6*Do2); rl = Sqrt[(6*Do2*po)/(omega*a)];

(*Functions needed*)
rn = Piecewise[{{0, ro <=  rl}, {ro*(0.5 - Cos[(ArcCos[1 - (2*rl*rl)/(ro^2)] - 2*Pi)/3]), ro  > rl}}];

p[r_] = Piecewise[{{po + con*(r^2 - ro^2 + 2*(rn^3)*(1/r - 1/ro)), r > rn} , {0 , r <= rn}}];

q[r_] = qm*((kme/(kme + p[r]))*(p[r]/(kmn + p[r]))  + (kmn/(kmn + p[r]))*j);

(*Equations defined*)
eqnDe = D[Ef1[r, t], t] - De*(D[Ef1[r, t], r, r] + (2/r)*(D[Ef1[r, t], r])) + q[r]*Ef1[r, t];

eqnBo = D[Eb1[r, t], t] - (Ef1[r, t])*q[r];

(*Solving equations*)
x = NDSolve[{eqnBo == 0, Eb1[r, 0] == 0, eqnDe == 0, Ef1[r, 0] == 0, Derivative[1, 0][Ef1][rn, t] == 0, Ef1[ro , t] == eo}, 
  Eb1, {r, rn, ro}, {t, 0, 14400}];

Then it's easy to plot Eb1 between rn and ro with this;

Plot[Eb1[r, 14400] /. x, {r, rn, ro}, PlotRange -> Automatic]

and I get this;

enter image description here

Yet when I try to flow Ef1 with a similar command I get absolutely nothing;

Plot[Ef1[r, 14400] /. x, {r, rn, ro}, PlotRange -> Automatic]

enter image description here

Similarly, I can't seem to evaluate Ef1 at a given value of $r$ or $t$ - any ideas what I'm doing wrong and why this is the case?

share|improve this question
add comment

1 Answer 1

up vote 1 down vote accepted

It seems that you're missing a solution of Ef1. Try

y = NDSolve[{eqnBo == 0, Eb1[r, 0] == 0, eqnDe == 0, Ef1[r, 0] == 0, 
Derivative[1, 0][Ef1][rn, t] == 0, Ef1[ro, t] == eo}, 
Ef1, {r, rn, ro}, {t, 0, 14400}];
Plot[Ef1[r, 14400] /. y, {r, rn, ro}, PlotRange -> Automatic]

enter image description here

share|improve this answer
1  
Thank you sir! Works absolutely perfectly. –  DRG Jul 7 at 15:49
    
@DRG No problem:-). –  Gregory Rut Jul 7 at 15:51
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.