# Linear ODEs with NDSolve

I am trying to solve some first order linear odes. My code is below:

eq={1. j1[z] +
10. u1[z] + (-0.04 - 0.02 λ[z]) w0[z] + 0.02 w1'[z] == 0,
(1.68333 + λ[z] (2. + 1. λ[z])) u0[z] + λ[z] ((2. - 0.5 λ[z]) u2[z] -
250. w1[z] - 25. ϕ0[z] + 25. ϕ1[z] + 0.5 u3'[z]) == 0,
100 u1[z] + ϕ1'[z] == 50 u3[z] + (2 + λ[z]) ϕ0[z],
(2 + λ[z]) j0[z] == 50 w1[z] + j1'[z],
u1[z] == u0'[z], u2[z] == u1'[z], u3[z] == u2'[z],
ϕ1[z] == ϕ0'[z], w1[z] == w0'[z], j1[z] == j0'[z]};

eqns = Rationalize[{eq, u0[0] == 0, u0[1] == 0,
u1[0]== 0, u1[1]== 0, w0[0] == 0, w0[1] == 0,
j1[0] == 0, j1[1] == 0, ϕ0[0] == 0,
ϕ0[1] == 0, λ'[z] == 0, u2[1] == 0}];

sol1 = NDSolve[
eqns, {u0, w0, λ}, {z, 0, 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {λ[0] == 20+ I},
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}
}(*, AccuracyGoal -> 10,PrecisionGoal -> 10*)
];

{λ[0] /. sol1}


There are 10 homogeneous equations and related 10 homogeneous boundary conditions. It is going to be an eigenvalue problem. For reasons unknown to me, the method doesn't converge. How can see where I'm going wrong while using the NDSolve? Are there any other ways to solve it. Please help. Thank you.

• I get some errors related to initial conditions. Have you tried other methods instead of the shooting method? Particularly, have you tried not specifying a method and letting NDSolve choose? Jun 15, 2022 at 16:26
• Consider also rewriting your equations completely to get rid of the Subscript expressions. They are a pretty-printing construct that sometimes cause misbehavior that is hard to track down. Jun 15, 2022 at 16:28
• @MarcoB Thank you. I have corrected the mistake. Yes, I meant eqns. I haven't tried it without giving a method. Jun 15, 2022 at 16:29
• @MarcoB I have made the changes. Jun 15, 2022 at 16:39
• Do you mean nonlinear instead of linear? (Technically the system is not homogenous. To be homogenous, all terms should have the same degree. At a glance, I see degree-1 and degree-2 terms. When a system is linear, the textbooks tend to emphasize the absence of a term without a dependent variable, which means all the terms present are degree-1.) Jun 16, 2022 at 0:21

This eigenvalue problem can be solved largely symbolically as follows. With the equations rationalized,

eq = {j1[z] + 10 u1[z] + w0[z] (-(1/25) - λ[z]/50) + w1'[z]/50 == 0,
u0[z] (168333/100000 + λ[z] (2 + λ[z])) + λ[z] (-250 w1[z] +
u2[z] (2 - λ[z]/2) - 25 ϕ0[z] + 25 ϕ1[z] + u3'[z]/2) == 0,
100 u1[z] + ϕ1'[z] == 50 u3[z] + (2 + λ[z]) ϕ0[z],
j0[z] (2 + λ[z]) == 50 w1[z] + j1[z],
u1[z] == u0'[z], u2[z] == u1'[z], u3[z] == u2'[z],
ϕ1[z] == ϕ0'[z], w1[z] == w0'[z], j1[z] == j0'[z]};


solve them symbolically with only the boundary conditions at z = 0 specified.

sol = DSolve[{eq /. λ[z] -> λ,
u0[0] == 0, u1[0] == 0, w0[0] == 0, j1[0] == 0, ϕ0[0] == 0},
{j0[z], j1[z], u0[z], u1[z], u2[z], u3[z], w0[z], w1[z], ϕ0[z], ϕ1[z]},
z] // Flatten


The solution is enormous and cannot be reproduced here. Nonetheless, the eigenvalues can be extracted by applying the z = 1 boundary conditions and constructing the array of coefficients of the five remaining constants of integration.

{u0[z], u1[z], w0[z], j1[z], ϕ0[z]} /. sol /. z -> 1;
co = Last@CoefficientArrays[%, {C[1], C[5], C[6], C[8], C[10]}] // Normal


Now, plot the determinant of this array as a function of λ to obtain estimates of the first several eigenvalues.

Plot[10^-15 Det[co], {λ, -500, 500}, PlotRange -> {-50, 50},
WorkingPrecision -> 30, ImageSize -> Large,
AxesLabel -> {λ, "Det"}, LabelStyle -> {15, Bold, Black}]


The first seven (in absolute value) are

N@FindRoot[10^-15 Det[co], {λ, -2.1, -1.9}, WorkingPrecision -> 30,
AccuracyGoal -> 3, PrecisionGoal -> 3, Evaluated -> False, Method -> "Brent"]
(* {λ -> -2.00116} *)


and, similarly,

(* {λ -> -15.3281} *)
(* {λ -> -35.2267} *)
(* {λ -> 110.078} *)
(* {λ -> 148.461} *)
(* {λ -> -201.539} *)
(* {λ -> -225.151} *)


The procedure just described works only for ODEs that can be solved using DSolve. The OP asked in a comment how to obtain eigenvalues for ODEs that can be solved only with NDSolve. The following works well.

The symbolic solution given above in effect seeks a linear combination of five independent solutions that satisfy the five boundary conditions at z = 1. This is possible only if the determinant of the coefficient matrix of those five solutions at z = 1 vanishes. That occurs only for isolated values of λ, determined using FindRoot. A similar process is feasible for numerical solutions. Obtain five linearly independent solutions satisfying the z = 0 boundary conditions by setting four of the other five dependent variables, {u2[z], u3[z], w1[z], j0[z], ϕ1[z]} to 0 there, and the other equal to 1. Then the corresponding values of {u0[z], u1[z], w0[z], j1[z], ϕ0[z]} comprise the rows of the coefficient matrix. For convenience, define

var = {u0[z], u1[z], w0[z], j1[z], ϕ0[z], u2[z], u3[z], w1[z], j0[z], ϕ1[z]};


Then the determinant of the coefficient matrix is given by

g[λt_] := Module[{v = Rationalize[λt, 0]}, Do[
bc0 = ReplacePart[ConstantArray[0, 10], n + 5 -> 1];
a[n] = NDSolveValue[{eq /. λ[z] -> v, Thread[(var /. z -> 0) == bc0]},
var[[;; 5]] /. z -> 1, {z, 0, 1}, WorkingPrecision -> 30], {n, 5}];
-10^-15 Det[{a[1], a[2], a[3], a[4], a[5]}]]


Plotting this function

Plot[g[λ], {λ, -500, 500}, PlotRange -> {-45, 45},
WorkingPrecision -> 30, MaxRecursion -> 1, ImageSize -> Large,
AxesLabel -> {λ, "Det"}, LabelStyle -> {15, Bold, Black}]


yields a figure indistinguishable to the eye from the plot given above, and the first seven eigenvalues are the same to five significant figures (except for the first eigenvalue which is the same only to three significant figures).

For completeness, corresponding eigenfunctions can be computed numerically without difficulty. With g defined as above but with WorkingPrecision -> 55 and MaxSteps -> 50000, compute the desired eigenvalue to greater accuracy, for instance,

tst = λ /. FindRoot[g[λ], {λ, -17, -15}, WorkingPrecision -> 55,
Evaluated -> False, Method -> "Brent"]
(* -15.328093316320908139126966813678620281397158355749 *)


Because the determinant of the coefficient matrix determined by g vanishes to high order, the coefficient matrix has a null space vector, which can be used as follows to determine the eigenfunctions.

NullSpace[Transpose[{a[1], a[2], a[3], a[4], a[5]}]] // Flatten;
bc0 = Join[ConstantArray[0, 5], Rationalize[%, 0]];
NDSolveValue[{eq /. λ[z] -> Rationalize[tst, 0], Thread[(var /. z -> 0) == bc0]},
var, {z, 0, 1}, WorkingPrecision -> 55];


Plots of the eigenfunction, first for the components with homogeneous boundary conditions and then for those with unspecified boundary conditions, are

Plot[Evaluate[%[[;; 5]]], {z, 0, 1}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {z}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[LineLegend[Head /@ var[[;; 5]],
LegendLayout -> "Row"], {.7, .1}]]
Plot[Evaluate[%%[[6 ;;]]], {z, 0, 1}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {z}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[LineLegend[Head /@ var[[6 ;;]],
LegendLayout -> "Row"], {.7, .1}]]


The accuracy with which the z = 1 boundary conditions are satisfied is

N[%%%[[;; 5]] /. z -> 1]
(* {1.06179*10^-10, 5.22517*10^-9, 5.09433*10^-8, -2.53322*10^-6, 2.59614*10^-7} *)


Incidentally, by inspection the smallest eigenvalue (in absolute value) is exactly λ = -2, and the corresponding eigenfunction is j0[z] a non-zero constant, and all other dependent variables identically zero.

Out of curiosity, I sought complex solutions and found two, plus their complex conjugates.

tst = λ /. FindRoot[g[λ], {λ, 300 + 100 I}, WorkingPrecision -> 55,
Evaluated -> False]
(* 238.2470380128362751886150009319777802016484715271720041 +
5.046719683312651069843507113048343649792611971592912802 I *)


The corresponding eigenfunctions, computed as before, are plotted using ReImPlot.

Another is

(* 440.9005698861846836113482607553640570450827531146712965 +
3.720425991419978280520680668723839826268342762230923495 I *)


Additional complex λ seem likely for larger Re[λ]. I attempted to find complex λ for Re[λ] < 100 but without success.

In a comment below, the OP asked assistance for modified boundary conditions, w1 instead of w0 vanishing at z = 0 asnd z = 1, and u1 changed to

m = 1.495; theta = -0.78;
bc = Rationalize[{u2[0] + m theta u1[0] == 0,
u2[1] + m theta/(1 + theta) u1[1] == 0}, 0]


For convenience, redefine var as

var = {u0[z], w1[z], j1[z], ϕ0[z], u1[z], u2[z], u3[z], w0[z], j0[z], ϕ1[z]};


Then, the definition of g is changed to incorporate the new boundary conditions.

g[λt_] := Module[{v = Rationalize[λt, 0]}, Do[
bc0 = ReplacePart[ConstantArray[0, 10], n + 5 -> 1];
bc0 = ReplacePart[bc0, 5 -> First[SolveValues[bc[[1]], u1[0]] /. u2[0] -> bc0[[6]]]];
a[n] = Flatten@NDSolveValue[{eq /. λ[z] -> v,
Thread[(var /. z -> 0) == bc0]}, {var[[;; 4]] /. z -> 1,
bc[[2, 1]]}, {z, 0, 1}, WorkingPrecision -> 30], {n, 5}];
-10^-19 Det[{a[1], a[2], a[3], a[4], a[5]}]]


and plotted.

Plot[g[λ], {λ, -500, 500}, PlotRange -> {-8, 6},
WorkingPrecision -> 30, MaxRecursion -> 3, ImageSize -> Large,
AxesLabel -> {λ, "Det"}, LabelStyle -> {15, Bold, Black}]


By inspection, the change in boundary conditions has no impact on the λ = -2 eigenvalue. Some others, determined as before, are

(* λ -> -17.0804038779156676304536707040 *)
(* λ -> -60.1390404619194112940862825698 *)
(* λ -> 45.4674271506238708342777922448 *)
(* λ -> 86.0968951335092628679561771189 *)


Computing these eigenvalues took as long as a minute each. As before, it is necessary to define

WorkingPrecision -> 55, MaxSteps -> 50000


in g to compute corresponding eigenfunctions. An example is

tst = λ /. FindRoot[g[λ], {λ, -18, -16},
WorkingPrecision -> 55, Evaluated -> False, Method -> "Brent"]

NullSpace[Transpose[{a[1], a[2], a[3], a[4], a[5]}]] // Flatten
bc0 = Join[ConstantArray[0, 5], Rationalize[%, 0]];
bc0 = ReplacePart[bc0,
5 -> First[SolveValues[bc[[1]], u1[0]] /. u2[0] -> bc0[[6]]]];
NDSolveValue[{eq /. λ[z] -> Rationalize[tst, 0],
Thread[(var /. z -> 0) == bc0]}, var, {z, 0, 1}, WorkingPrecision -> 55];

Plot[Evaluate[%[[;; 5]]], {z, 0, 1}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {z}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[LineLegend[Head /@ var[[;; 5]],
LegendLayout -> "Row"], {.7, .1}]]
Plot[Evaluate[%%[[6 ;;]]], {z, 0, 1}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {z}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[LineLegend[Head /@ var[[6 ;;]],
LegendLayout -> "Row"], {.7, .1}]]
Join[N[%%%[[;; 4]] /. z -> 1],
{N[bc[[2, 1]] /. ({u1[1] -> %%%[[5]], u2[1] -> %%%[[6]]}) /. z -> 1]}]


(* {-4.74378*10^-11, -1.11753*10^-6, 1.12926*10^-6, -1.15817*10^-7, -1.02345*10^-7} *)


The final line shows that the new boundary conditions are well satisfied.

• Thank you for your wonderful method. It works fine. However, I have a query. Can it be used when the coefficients are a function of the variable (in my case z)? I have a problem with those types of coefficients as well. The problem with NDSolve is it requires an extra condition on the solution to determine the eigenvalue. Jun 16, 2022 at 7:50
• @bbgodfrey, please notice the typo: The OP has w1[z] == w0'[z] . Your answer uses w1[z] == w0[z] Jun 16, 2022 at 14:33
• @GAURAVMAURYA As requested, I have added a solution using NDSolve instead of DSolve. Jun 16, 2022 at 23:52
• @GAURAVMAURYA And, now the corresponding eigenfunctions. Jun 18, 2022 at 20:11
• @MichaelE2 The symbolic solution is quite slow, which I partially circumvented with AccuracyGoal -> 3, PrecisionGoal -> 3. The numerical solution is both faster and more accurate. The computation in the second addendum takes about a minute. This is true both for the specific eigenvalue shown and for several others. Jun 19, 2022 at 14:55

## Update

As mentioned by bbgodfrey and by Michael E2, the code I had proposed below only returns a trivial constant solution, as the following indicates:

Through[sol1[[1, All, 2]]["ValuesOnGrid"]] // Chop

(* Out:
{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{20.00000000000000000000000 + 1.000000000000000000000000 I,
20.00000000000000000000000 + 1.000000000000000000000000 I,
[...]
20.00000000000000000000000 + 1.000000000000000000000000 I,
20.00000000000000000000000 + 1.00000000000000000000000 I}}
*)


Replacing your original Subscript expressions helped a lot. Once those were removed, your original code produced an error from FindRoot, suggesting that a higher working precision could be needed. After some trial and error, a WorkingPrecision -> 25 setting seems to work.

In the following code I also replaced your machine-precision constants (e.g. 0.04) with arbitrary-precision ones (e.g. 4/100), and combined equations and boundary conditions for clarity.

eqandbc = {
j1[z] + 10 u1[z] + (-4/100 - 2/100 λ[z]) w0[z] + 2/100 w1'[z] == 0,
(101/60 + λ[z] (2 + λ[z])) u0[z] + λ[z] ((2 - 1/2 λ[z]) u2[z] -
250 w1[z] - 25 ϕ0[z] + 25 ϕ1[z] + 1/2 u3'[z]) == 0,
100 u1[z] + ϕ1'[z] == 50 u3[z] + (2 + λ[z]) ϕ0[z],
(2 + λ[z]) j0[z] == 50 w1[z] + j1'[z],
u1[z] == u0'[z], u2[z] == u1'[z], u3[z] == u2'[z],
ϕ1[z] == ϕ0'[z], w1[z] == w0'[z], j1[z] == j0'[z],

u0[0] == 0, u0[1] == 0,
u1[0] == 0, u1[1] == 0,
w0[0] == 0, w0[1] == 0,
j1[0] == 0, j1[1] == 0,
ϕ0[0] == 0, ϕ0[1] == 0,
λ'[z] == 0, u2[1] == 0
};

sol1 =
NDSolve[
eqandbc, {u0, w0, λ}, {z, 0, 1},
Method -> {
"Shooting",
"StartingInitialConditions" -> {λ[0] == 20 + I},
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}
},
WorkingPrecision -> 25
]


• Thank you so much. So, basically I need to set working precision to get it working and it will work for other bcs as well. Jun 16, 2022 at 1:36
• @GAURAVMAURYA That's likely, although performance may sometimes depend on the exact boundary conditions. But just try and see. Jun 16, 2022 at 1:37
• I noticed that the value of λ returned by your computation was essentially identical to your initial guess, λ[0] == 20 + I. So I tried other initial guesses, and the code again returned the initial guess, which seems suspicious. Jun 16, 2022 at 4:51
• @MarcoB Thank you. Jun 16, 2022 at 5:02
• As @bbgodfrey hinted, this this does not appear to be an effective solution. The solution produce is in effect the trivial (zero) solution, which would be apparent if the answer showed the plots. Zero is a solution for every proposed eigenvalue λ, but it is not an eigenfunction or a solution to the eigenvalue problem. Jun 19, 2022 at 14:49

I have a package for solving eigenvalue boundary value problems using the Compound Matrix Method with the Evans function, which I'll use here. The package is available on my GitHub (which has a notebook with examples), more details are in my other answers to questions here, and a good introduction to the method is in this pdf.

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod", "Site" ->
"http://raw.githubusercontent.com/paclets/Repository/master"]


Basically the Compound Matrix Method takes an $$n$$ by $$n$$ eigenvalue problem of the form $$\mathbf{y}' = A(x, \lambda) \cdot \mathbf{y}, \quad a \leq x \leq b, \\ B(x,\lambda) \mathbf{y} = \mathbf{0}, \quad x=a, \\ C(x,\lambda) \mathbf{y} = \mathbf{0}, \quad x=b,$$ and converts it to a larger system of determinants that satisfy a different matrix equation $$\mathbf{\phi}' = Q(x, \lambda) \mathbf{\phi}.$$ This removes a lot of the stiffness from the equations, as well as being able to also remove the exponential growth terms that dominate away from an eigenvalue.

This is quite similar to bbgodfrey's solution of finding the determinant, but is more numerically stable.

I've taken your equations back to a natural form, without the first order substitutions, including just using $$u, w, j, \phi$$. The following function will sort all that out for you.

eqns = {j'[z] + 10 u'[z] + (-4/100 - 2/100 λ) w[z] + 2/100 w''[z] == 0,
(101/60 + λ (2 + λ)) u[z] + λ ((2 - 1/2 λ) u''[z] - 250 w'[z] -
25 ϕ[z] + 25 ϕ'[z] + 1/2 u''''[z]) == 0,
100 u'[z] + ϕ''[z] == 50 u'''[z] + (2 + λ) ϕ[z],
(2 + λ) j[z] == 50 w'[z] + j''[z]};

bcs = {u[0] == 0, u[1] == 0, u'[0] == 0, u'[1] == 0, w'[0] == 0,
w'[1] == 0, j'[0] == 0, j'[1] == 0, ϕ[0] == 0, ϕ[1] == 0};


We write the system in the matrix form $$\mathbf{y}' = \mathbf{A} \cdot \mathbf{y}$$, with the function ToMatrixSystem:

sys = ToMatrixSystem[eqns, bcs, {u, j, w, ϕ}, {z, 0, 1}, λ]


Then we use the "Evans Function", an analytic function whose roots correspond to the eigenvalues of the original system, reducing the problem to finding roots of a smooth function of one variable.

We can evaluate this function at a point, for instance $$\lambda = -3$$:

Evans[-3, sys]
(* -4.41549*10^-10 *)


This value is not zero, so $$\lambda = -3$$ is not an eigenvalue.

We can plot this function, and see that it is smooth, with a (double) root around 2:

Plot[Evans[\[Lambda], sys], {\[Lambda], -5, 0}


And we can use FindRoot to find a few roots:

FindRoot[Evans[λ, sys], {λ, -3}]
(* {λ -> -2.} *)
FindRoot[Evans[λ, sys], {λ, -15}]
(* {λ -> -14.5244} *)
FindRoot[Evans[λ, sys], {λ, -100}]
(* {λ -> -134.338} *)


(This isn't exhaustive by any means, just a few that I found immediately).

I'm not sure why these is different to bbgodfrey's values though, I haven't dug into that. Perhaps I have made a mistake with the equations.

This method should work with coefficients as functions of $$z$$. I know it isn't Mathematica, but I would also say that I found the MATLAB plugin chebfun to be excellent for this kind of problem. As requested, a scrap of chebfun MATLAB code for a 10th order system (my MATLAB skill is very limited, I've taken some old code and removed constants but can't test this still works).

% Assign the equation to two chebops N and B such that N(u) = lambda*B (u).
dom = [0 1];
N = chebop(@(x,u1,u10,u2,u3,u4,u5,u6,u7,u8,u9) [diff(u1)-u2-u7-u9;diff(u2)-u3;diff(u3)-u4;diff(u4)-u2-u3-u6-u7-u9;diff(u5)-u6;diff(u6)-u4-u5;diff(u7)-u8;diff(u8)-u2-u7-u8-u9;diff(u9)-u10;diff(u10)-u2-u7-u9-u10], dom);
B = chebop(@(x,u1,u10,u2,u3,u4,u5,u6,u7,u8,u9) [0.*u1;0.*u2;0.*u3;u1;0.*u5;0.*u6;0.*u7;1.*u7;0.*u9;1.*u9], dom);

% Assign boundary conditions to the chebop.
N.bc = @(x,u1,u10,u2,u3,u4,u5,u6,u7,u8,u9) [u2(0); u3(0); u5(0); u8(0); u10(0); u2(1); u3(1); u5(1); u8(1); u10(1)];

%% Setup preferences for solving the problem.
% Create a CHEBOPPREF object for passing preferences.
% (See 'help cheboppref' for more possible options.)
options = cheboppref();
options.discretization = 'values';

%% Solve the eigenvalue problem.
[V, D] = eigs(N, B, k, options,'SM');
D = diag(D);

• Thank you for this function. I have a query. What if there resultant values are complex numbers? I mean, you have only plotted the real part. Jun 28, 2022 at 12:00
• The Evans function is in general a function that maps from the complex plane to the complex plane. If you have entirely real values in your system (including boundary conditions) and give a real eigenvalue I think that the Evans function will be real if I remember correctly. But you can give e.g. Evans[1+I,sys] and it will give you a complex value, and FindRoot can find complex values. Jun 28, 2022 at 12:23
• I would suggest taking a look at the examples that I include in the github, i.e. file CMMExamples.nb from github.com/SPPearce/CompoundMatrixMethod . I have some examples of complex valued systems and eigenvalues there. Note that you can have complex eigenvalues even given an entirely real system, I just think that if you try a real value of $\lambda$ then $E(\lambda)$ will be real if the system is real. Jun 28, 2022 at 12:24
• thank you so much. The method works fine to find the eigenvalues. But when I try to use the found eigenvalues in NDSolve to find the eigenvectors, it doesn't work. If you know the reason, please share. Thank you once again. Jun 30, 2022 at 9:14
• I have added some chebfun` code that may be useful to you, although I can't test it. My package doesn't easily give you the eigenvectors I'm afraid. It is tricky to set that up for the general case, so I have never got round to coding that up (if it is indeed possible). Jun 30, 2022 at 10:28