From a system of PDEs where i used the following ansatz: $$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$$. $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$ So, $$\theta_w(x,y) = e^{-\beta_h x} F'(x) e^{-\beta_c y} G'(y)$$ I have two linear ODEs of which one i am mentioning here $$ \lambda_hF'''-2\lambda_h\beta_hF''+((\lambda_h\beta_h-1)\beta_h-\mu)F'+\beta_h^2F=0 $$ with the b.c. $$ F(0)=0 $$ $$ \frac{F''(0)}{F'(0)}=\beta_h $$

$$ \frac{F''(1)}{F'(1)}=\beta_h$$

Here $\mu$ is the separation constant

A general solution of this ODE will be of the form:

$$ F(x) = \sum_k C_k e^{-\delta_k(\mu)x}$$ $k$ attains values $1,2,3$ because of the cubic nature of the characterstic equation.

I need to find the eigenvalues of this ODE. I know i need to consider $\mu><=$$0$. A second order ODE has a well defined form of its roots which is not the case here.

I will mention here that i am not too well accustomed to Mathematica But is there any way through which these eigen values could be determined using Mathematica ? I am at a dead end here.

(Finally i need to repeat the complete procedure for $G$ and then proceed to find $\theta_w$)

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Jan 14, 2019 at 17:07
  • 4
    $\begingroup$ People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful $\endgroup$
    – Michael E2
    Jan 14, 2019 at 17:07
  • 1
    $\begingroup$ Also the boundary condition at 1 is identical on the top and the bottom of the fraction $\endgroup$
    – SPPearce
    Jan 14, 2019 at 18:38
  • $\begingroup$ I can solve it numerically for a given value of $\beta$ and $\lambda$ (assuming it should be $F''(1)/F'(1)=\beta$). Depending on their values, both positive and negative eigenvalues can exist. $\endgroup$
    – SPPearce
    Jan 14, 2019 at 19:46
  • 1
    $\begingroup$ Try solving the ODE and two of the three boundary conditions with DSolve. Then apply the third boundary condition to the resulting output to obtain an algebraic equation for the eigenvalues. $\endgroup$
    – bbgodfrey
    Jan 14, 2019 at 20:08

2 Answers 2


I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions or the github for some more details.

First we install the package (only need to do this the first time):

    "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Then we first need to turn the ODEs into a matrix form $\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$, using my function ToMatrixSystem:

sys = ToMatrixSystem[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0}, 
         {F[0] == 0, F''[0] == β F'[0], F''[1] == β F'[1]}, F, {x, 0, 1}, μ]

The object sys contains the matrix $\mathbf{A}$, as well as similar matrices for the boundary conditions and the range of integration.

Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $\mu$ provided that $\lambda$ and $\beta$ are specified; this is an analytic function whose roots coincide with eigenvalues of the original equation.

Evans[1, sys /. {λ -> 1, β -> 2}]
(* -0.023299 *)

We can also plot this:

Plot[Evans[μ, sys /. {λ -> 1, β -> 1}], {μ, -300, 300}, AxesLabel -> {"μ", "E(μ)"}]

enter image description here

For these parameter values you can see there are six eigenvalues within the plot region, the function findAllRoots (by user Jens, available from this post) will give you them all:

findAllRoots[Evans[μ, sys /. {λ -> 1, β -> 1}], {μ, -300, 300}]
(* {-247.736, -158.907, -89.8156, -40.4541, -10.7982, -0.635232} *)

For the values $\lambda=0.02, \beta=10$, it helps to remove the normalisation that is applied by default in Evans, and I've also plotted the values given by bbgodfrey's solution:

Plot[Evans[μ, sys /. {λ -> 1/50, β -> 10}, NormalizationConstants -> 0], 
    {μ, -21, -2}, Epilog -> 
  Point[{#, 0} & /@ {-20.88, -17.48, -14.36, -11.53, -9.03, -6.92, -5.26, -4.08, -3.16}]]

enter image description here

It looks like there is an infinite set of solutions for negative $\mu$ for this set of parameters. For $\beta=0$ the eigenvalues are $-\lambda (n \pi)^2$, the general behaviour as $\mu->-\infty$ looks the same.

As you vary $\beta$ and $\lambda$ you can see how the function changes and eigenvalues can move and coalesce.

Note there is a discrepancy with bbgodfrey's results, which missed the root at $-3.392$ and gave a spurious one at $-3.165$. That isn't to say that you can't manage to get spurious solutions by applying FindRoot to this technique if you aren't careful.

Additionally, as the Evans function is (complex) analytic, you can use the Winding Number to find how many roots are within a contour. I have some functions for this, so you can see that there are only these 9 roots within the circular contour of size 21 located at the origin:

PlotEvansCircle[sys /. {λ -> 1/50, β -> 10}, ContourRadius -> 21, 
   nPoints -> 500, Joined -> True]

enter image description here

Here the left hand side is a circle in the complex plane, and the right hand side is the Evans function applied to the points on that circle. The winding number is how many times the right hand contour goes round the origin. There is also a PlotEvansSemiCircle for seeing if there is a root in the positive half-plane (as that is a common question for instability).

To get the eigenfunctions, replace the right BC with $F'(0)$ to some arbitrary value and use NDSolve. This is at high precision to check the boundary conditions are being met:

Block[{β = 10, λ = 1/50},
 μc = μ /. FindRoot[Evans[μ, sys, NormalizationConstants -> 0, 
      WorkingPrecision -> 50], {μ, -3}, WorkingPrecision -> 50, 
     AccuracyGoal -> 15, PrecisionGoal -> 15] // Quiet;
 sol = First[NDSolve[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0, 
   F[0] == 0, F''[0] == β F'[0], F'[0] == 1/10} /. μ -> μc, F, 
   {x, 0, 1}, WorkingPrecision -> 50, AccuracyGoal -> 30, PrecisionGoal -> 30]];
 Plot[(F[x]/F[1]) /. sol, {x, 0, 1}, 
  PlotLabel -> 
   "Eigenfunction for μ=" <> ToString@Round[μc, 0.0001] <> ", with BC error: " <> 
    ToString@Round[((F''[1] - β F'[1]) /. sol), 0.0001], 
  AxesLabel -> {"x", "F(x) (scaled)"}, PlotRange -> All]


The eigenfunctions get increasingly oscillatory as the eigenvalue gets more negative:

roots = Reverse@findAllRoots[Evans[μ, sys /. {λ -> 1/50, β -> 10}, 
    NormalizationConstants -> 0], {μ, -200, -2}]
Plot[Evaluate[Table[(F[x]/F[1]) /. 
    NDSolve[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0, 
 F[0] == 0, F''[0] == β F'[0], F'[0] == 1/1000} /. {λ -> 1/50, β -> 10}, 
     F, {x, 0, 1}], {μ, roots}]], {x, 0, 1}, PlotRange -> All, 
 AxesLabel -> {"x", "F(x), scaled"}]

enter image description here

The method will work for higher order systems (up to 10th order generally), and does not require DSolve to be able to solve the underlying ODE.

  • $\begingroup$ I should probably note that I don't currently have support for getting the eigenfunctions out of my package, which it looks like you'll probably need to proceed. $\endgroup$
    – SPPearce
    Jan 15, 2019 at 10:13
  • $\begingroup$ This was great and very useful. My only doubt was regarding the range of $\mu$ selected, which i guess i need to choose judiciously. Further since the general solution of $F$ is of the form $ F(x) = \sum_k C_k e^{-\delta_k(\mu)x}$, the eigen functions should be of the form $e^{-\delta_k(\mu_n)t}$ The $\delta_k$ is the root of the characteristic equation $\endgroup$
    – Avrana
    Jan 15, 2019 at 11:56
  • $\begingroup$ @IndrasisMitra, added how to get the eigenfunctions and some examples. I would say not to accept an answer immediately, and wait to see if any other answers are forthcoming. $\endgroup$
    – SPPearce
    Jan 15, 2019 at 13:28
  • 1
    $\begingroup$ can the increasing oscillatory nature of eigenfunctions as the eigenvalues get more negative be inferred in a way that $F$ does not converge ? $\endgroup$
    – Avrana
    Jan 15, 2019 at 18:20
  • $\begingroup$ also i see that you have numerically solved the ODE with an arbitrary value of $F'(0)$ to find the eigenfuctions. I fail to grasp the logic behind this step $\endgroup$
    – Avrana
    Jan 15, 2019 at 18:43

Error corrected by deleting Numerator from definition of s.

This problem can be transformed symbolically to an algebraic equation by the method described in my earlier comment.

eq = ll f'''[x] - 2 ll b f''[x] + ((ll b - 1) b - m) f'[x] + b^2 f[x] == 0;
s = Simplify[DSolveValue[{eq, f''[0] == b f'[0], f''[1] == b f'[1]}, f[x], x]
    /. x -> 0 /. C[1] -> 1] 

The resulting expression in terms of Root is too lengthy to be reproduced here and must, in any case, be solved numerically to provide values for m. For the parameters given by the OP in a comment,

Table[FindRoot[s /. {ll -> 1/50, b -> 10}, {m, m0}] // Chop // Values, 
    {m0, -21, -1, 1/10}] // Flatten;
Union[%, SameTest -> (Abs[#1 - #2] < 10^-4 &)]

(* {-20.8885, -17.4846, -14.3644, -11.5367, -9.03846, -6.92843, 
    -5.26567, -4.08628, -3.39259} *)

There appear to be no larger values of m for this set of parameters. There are, of course, many more-negative values, perhaps an infinite number.

Here are a few more cases. For {ll -> 1/10, b -> 10},

(* {-13.9638, -7.91839, -3.94417, -1.94864} *)

and for {ll -> 1/1000, b -> 10},

(* {-20.5949, -19.94, -19.3031, -18.684, -18.0821, -17.4972, -16.9288, 
    -16.3763, -15.8388, -15.3158, -14.806, -14.3085, -13.822, -13.3451, 
    -12.8764, -12.4145, -11.9582, -11.5069, -11.0604, -10.6193, -10.1852, 
    -9.76044, -9.34825, -8.95222, -8.57622, -8.22411, -7.89954, -7.6058, 
    -7.34579, -7.12193, -6.93623, -6.79024, -6.68516, -6.62179} *)
  • $\begingroup$ Thanks a lot for this help. Now i have some idea. Some points i need to clarify. For the case $\lambda_h=0.1$ and $\beta_h=10$, you have posted a set of eigenvalues that come out. You also mention that there could be infinitely many more negative values. Any advice on how many of these eigen values should one consider to generate a solution ? I know the query is not related to mathematica here. Thanks again. $\endgroup$
    – Avrana
    Jan 15, 2019 at 3:48
  • $\begingroup$ @IndrasisMitra There is no good way to know without seeing the actual sum you are constructing. Is it a sum of eigenfunctions for the ODE in your question, or something else? In general, if the function to be approximated by a sum is reasonable uniform, then probably only a half-dozen or so eigenvalues are needed. $\endgroup$
    – bbgodfrey
    Jan 15, 2019 at 5:14
  • $\begingroup$ I have edited my original post to show the type of series i expect. Yes, it is a sum of eigenfunctions from the ODE i posted. $\endgroup$
    – Avrana
    Jan 15, 2019 at 5:29
  • $\begingroup$ I'm getting a catastrophic precision loss for FindRoot at your largest root for {ll -> 1/50, b -> 10}: s /. {ll -> 1/50, b -> 10} /. FindRoot[s /. {ll -> 1/50, b -> 10}, {m, -3}, WorkingPrecision -> 20] gives me no significant digits to display (I don't believe it to be a root). $\endgroup$
    – SPPearce
    Jan 15, 2019 at 9:33
  • 1
    $\begingroup$ @KraZug I now see that it is not a root. My error was to use only the Numerator, which I plan to correct soon. Thanks. $\endgroup$
    – bbgodfrey
    Jan 15, 2019 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.