Solving boundary value problem with coupled odes at interface

I am trying to get the eigenvalues of the following differential system

eq1 = -2*\[Psi]'[r]/r^3 - k^2*\[Psi]''[r] + 2*\[Psi]''[r]/r^2 - k^2*(-k^2*\[Psi][r] - \[Psi]'[r]/r + \[Psi]''[r]) - \[Psi]'''[r]/r - (-k^2*\[Psi]'[r] + \[Psi]'[r]/r^2 - \[Psi]''[r]/r + \[Psi]'''[r])/r + \[Psi]''''[r] == 0;
eq2 = -k^2*\[Phi]A[r] + \[Phi]A'[r]/r + \[Phi]A''[r] == 0;
eq3 = -k^2*\[Phi]F[r] + \[Phi]F'[r]/r + \[Phi]F''[r] == 0;

with some functions used in the system

bw[r_] = -(1/4)*(r^2 - 1/4 - 2*Log[2*r]); bq = (\[Epsilon] - 1)/Log; b\[Phi]A[r_] = Log[2*r]/Log; b\[Phi]F[r_] = Log[2*r]/Log;

The odes are subjected to the boundary conditions at rL = 1/2 and rR = 2:

\[Psi][rL] == 0; \[Psi]'[rL] == 0; \[Phi]F[rL] == 0;
bcR = \[Phi]A[rR] == 0;

as well as the matching condtion at r=1:

mbc1 = k^2*\[Psi]*(c - bw) + \[Psi]''*(c - bw) - \[Psi]'*(c - bw) - \[Psi] == -I*k*2*bq*(\[Psi]*b\[Phi]A' + \[Phi]A*(c - bw));
mbc2 = k*\[Psi]*(1 - k^2) - 2*I*k^2*\[Psi]*(c - bw) + I*(3*k^2 - 1)*\[Psi]'*(c - bw) - I*\[Psi]'''*(c - bw) + I*\[Psi]''*(c - bw) == -2*k*(b\[Phi]A'*(b\[Phi]A''*\[Psi] + \[Phi]A'*(c - bw)) - \[Epsilon]*b\[Phi]F'*(b\[Phi]F''*\[Psi] + \[Phi]F'*(c - bw)));
mbc3 = (-I*k*c + I*bw*k)*(\[Epsilon]*(b\[Phi]F''*\[Psi] + \[Phi]F'*(c - bw)) - (b\[Phi]A''*\[Psi] + \[Phi]A'*(c - bw))) + I*k*(-\[Psi] + \[Psi]')*(c - bw)*bq == 5*(b\[Phi]A''*\[Psi] + \[Phi]A'*(c - bw)) - 5*(b\[Phi]F''*\[Psi] + \[Phi]F'*(c - bw));
mbc4 = \[Phi]F*(c - bw) + b\[Phi]F'*\[Psi] == \[Phi]A*(c - bw) + b\[Phi]A'*\[Psi];

in which c is a complex eigenvalue in general, k and \[Epsilon] are parameters.

Note:

1. the function bq and the b.c.s mbc2 and mbc3 include the parameter \[Epsilon];

2. eq1 and eq3 are defined in rL<=r<=1 while eq2 is defined in 1<=r<=rR;

3. eq2 and eq3 have the same form and both have a general (explicit) solution: \[Phi][r] == C1*BesselI[0, k*r] + C2*BesselK[0, k*r]. However, in this problem, I'd like to solve them numerically.

My aim is to calculate the eigenvalue c for a set of \[Epsilon] and k. I have tried to use the package developed by @SPPearce since it can deal with a similar problem with an interface. The main difference is that in my problem the eigenvalue only appears in the b.c.s, while in that problem the eigenvalue appears in both odes. I have also noted the package can cope with eigenvalue dependent b.c.s, which can be invoked as follows:

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

sys[k_, \[Epsilon]_] = With[{k = k, \[Epsilon] = \[Epsilon]}, ToMatrixSystem[{eq1, eq2, eq3}, {\[Psi][rL] == 0, \[Psi]'[rL] == 0, \[Phi]F[rL] == 0, bcR, mbc1, mbc2, mbc3, mbc4}, {\[Psi], \[Phi]F, \[Phi]A}, {r, rL, 1, rR}, c]]

Note that according to the problem the independent variable should be given as {r, rL, 1, rR}, however, ToMatrixSystem returns the system unevaluated and Plot[Evans[c, sys[1, 5]], {c, 1, 3}] gives null. With independent variable specified as {r, rL, 1} instead, it seems to put the equations into matrix form as required by this method. But Plot[Evans[c, sys[1, 5]], {c, 1, 3}] gives many errors. I understood that the problem could be converted into a solvability condition for a matrix problem: M x=0, in which we may require Det[M]==0 for non-trivial solutions. But in general for odes without explicit solutions, I'd like to solve the problem numerically.

This package works well for a problem with two coupled odes and for the above-mentioned problem with an interface, both of which look more complicated than mine. I don't understand why it does not work. I would be very thankful if anybody could suggest how to solve this problem.

Rather than trying to hack the package referenced in the question, I found it easier to solve the problem from first principles. With quantities as defined in the question, the solution is

csolve[k0_, ϵ0_] := Module[{sψ1, sψ2, sA, sF, rules, a1, a2, a3, a4},
sψ1 = NDSolveValue[{eq1 /. k -> k0, ψ[rL] == 0, ψ'[rL] == 0, ψ''[rL] == 1,
ψ'''[rL] == 0}, {ψ, ψ', ψ'', ψ'''}, {r, rL, 1}];
sψ2 = NDSolveValue[{eq1 /. k -> k0, ψ[rL] == 0, ψ'[rL] == 0, ψ''[rL] == 0,
ψ'''[rL] == 1}, {ψ, ψ', ψ'', ψ'''}, {r, rL, 1}];
sF = NDSolveValue[{eq3 /. k -> k0, ϕF[rL] == 0, ϕF'[rL] == 1},
{ϕF, ϕF'}, {r, rL, 1}];
sA = NDSolveValue[{eq2 /. k -> k0, ϕA[rR] == 0, ϕA'[rR] == 1},
{ϕA, ϕA'}, {r, 1, rR}];
rules = Join[Array[(D[ψ[r], {r, # - 1}] -> a1 sψ1[[#]] + a2 sψ2[[#]]) &, 4],
Array[(D[ϕF[r], {r, # - 1}] -> a3 sF[[#]]) &, 2],
Array[(D[ϕA[r], {r, # - 1}] -> a4 sA[[#]]) &, 2]] /. r -> 1;
NSolveValues[CoefficientArrays[Subtract @@@ {mbc1, mbc2, mbc3, mbc4} /. rules
/. k -> k0 /. ϵ -> ϵ0, {a1, a2, a3, a4}] // Last // Det, c]]

A typical solution is

csolve[4., .1]
(* {0.158348 - 2.1338 I, 0.159074, 0.159074, 0.159074, 0.169511 - 0.275478 I} *)

In general, two of the solutions are distinct and vary with k and ϵ. The other three are unchanging, equal to 0.159074.

The solution takes advantage of the fact that the equations, along with their boundary conditions, are linear and homogeneous. Hence, two independent solutions for eq1 satisfying the boundary condition at rL can be computed using NDSolveValue with ψ''[rR] and ψ'''[rL] specified arbitrarily. Note that only {ψ, ψ', ψ'', ψ'''} are specified as output, because only they are needed to compute c. Values for {ϕF, ϕF'} and {ϕA, ϕA'} are determined similarly. These quantities then are substituted into {mbc1, mbc2, mbc3, mbc4} to obtain four equations for the four coefficients of the ODE solutions previously computed. Because these linear equations are homogeneous, they have nontrivial solutions only if the determinant of the terms of the equations vanish, which determines c.

This procedure can be generaized without difficulty to other systems of linear, homogeneous ODEs.

• Thank you very much! Which version are you using? Why did you apply Last in CoefficientArrays to exclude nonhomogeneous constants in the coefficient matrix? Dec 27 '21 at 3:11
• @Nobody, there are no inhomogeneous constants in the coefficient matrix, as can be seen from CoefficientArrays[Subtract @@@ {mbc1, mbc2, mbc3, mbc4} /. rules /. k -> 1 /. ϵ -> 1, {a1, a2, a3, a4}] // First // Normal. I am using Mathematica 13.0 with Windows 10. Dec 27 '21 at 4:15
• Some observations: 1. in general the system has 5 solutions, 3 of which are the same fixed real number and the other 2 are distinct complex numbers varying with k and eps; 2. sometimes, it has 4 solutions (e.g. for k=0, eps=1), three of which are the fixed real number, the other is complex varying with k and eps. The fixed real and varying complex solutions do not depend on the arbitrary bcs. It suggests than the fixed real solution may be fictitious. If that is the case, should it be deleted, if not, might it represent an intrinsic characteristic of the system? Dec 27 '21 at 13:01
• In answer to your two questions, the results do not depend on the arbitrary b.c.s so long as the resulting NDSolve solutions are linearly independent, and they are. To verify the code in my answer, I copied it from my answer and compared the results with those of the symbolic solution. They agree within the precision of the computation. The fixed three identical values of c also occur in the symbolic solution, so they appear to be correct. By the way, increasing the precision of the numerical results for c, which is only about 10^-6, is not difficult. Dec 27 '21 at 21:14
• The three-fold constant root is given by c -> 1/16 (-3 + 8 Log). In other words, c -> bw, which may be spurious. Dec 28 '21 at 2:43