Solving boundary value problem with coupled odes at interface

Question

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[4]; b\[Phi]A[r_] = Log[2*r]/Log[4]; b\[Phi]F[r_] = Log[2*r]/Log[4];

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][1]*(c - bw[1]) + \[Psi]''[1]*(c - bw[1]) - \[Psi]'[1]*(c - bw[1]) - \[Psi][1] == -I*k*2*bq*(\[Psi][1]*b\[Phi]A'[1] + \[Phi]A[1]*(c - bw[1]));
mbc2 = k*\[Psi][1]*(1 - k^2) - 2*I*k^2*\[Psi][1]*(c - bw[1]) + I*(3*k^2 - 1)*\[Psi]'[1]*(c - bw[1]) - I*\[Psi]'''[1]*(c - bw[1]) + I*\[Psi]''[1]*(c - bw[1]) == -2*k*(b\[Phi]A'[1]*(b\[Phi]A''[1]*\[Psi][1] + \[Phi]A'[1]*(c - bw[1])) - \[Epsilon]*b\[Phi]F'[1]*(b\[Phi]F''[1]*\[Psi][1] + \[Phi]F'[1]*(c - bw[1])));
mbc3 = (-I*k*c + I*bw[1]*k)*(\[Epsilon]*(b\[Phi]F''[1]*\[Psi][1] + \[Phi]F'[1]*(c - bw[1])) - (b\[Phi]A''[1]*\[Psi][1] + \[Phi]A'[1]*(c - bw[1]))) + I*k*(-\[Psi][1] + \[Psi]'[1])*(c - bw[1])*bq == 5*(b\[Phi]A''[1]*\[Psi][1] + \[Phi]A'[1]*(c - bw[1])) - 5*(b\[Phi]F''[1]*\[Psi][1] + \[Phi]F'[1]*(c - bw[1]));
mbc4 = \[Phi]F[1]*(c - bw[1]) + b\[Phi]F'[1]*\[Psi][1] == \[Phi]A[1]*(c - bw[1]) + b\[Phi]A'[1]*\[Psi][1];

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.

Answer 1

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}, {ψ[1], ψ'[1], ψ''[1], ψ'''[1]}, {r, rL, 1}]; 
  sψ2 = NDSolveValue[{eq1 /. k -> k0, ψ[rL] == 0, ψ'[rL] == 0, ψ''[rL] == 0, 
    ψ'''[rL] == 1}, {ψ[1], ψ'[1], ψ''[1], ψ'''[1]}, {r, rL, 1}]; 
  sF = NDSolveValue[{eq3 /. k -> k0, ϕF[rL] == 0, ϕF'[rL] == 1}, 
    {ϕF[1], ϕF'[1]}, {r, rL, 1}]; 
  sA = NDSolveValue[{eq2 /. k -> k0, ϕA[rR] == 0, ϕA'[rR] == 1}, 
    {ϕA[1], ϕA'[1]}, {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 {ψ[1], ψ'[1], ψ''[1], ψ'''[1]} are specified as output, because only they are needed to compute c. Values for {ϕF[1], ϕF'[1]} and {ϕA[1], ϕA'[1]} 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.