I must say that I don't much understand the motivation coming from number theory, but your question about rational solutions of ODEs has a definite answer, provided the equation has polynomial or rational coefficients. A standard reference is
Abramov, S. A., Rational solutions of linear differential and difference equations with polynomial coefficients, U.S.S.R. Comput. Math. Math. Phys. 29, No. 6, 7-12 (1989). ZBL0719.65063.
The main idea is straight forward and I will summarize it here. The solution of a homogeneous linear ODE with rational coefficients (which can always be made polynomial by clearing denominators) can have a pole (including the poles at infinity) only at a regular singular point, and only if the Frobenius method reveals integer negative critical exponents. Taking the most negative integer critical exponents at each regular singular point gives you a universal polynomial denominator $Q$, so that any rational solution can be written as $P/Q$, with $P$ polynomial. The Frobenius analysis at infinity also gives you an upper bound on the degree of $P$, which makes $P/Q$ an ansatz for the most general rational solution with finitely many parameters. Now, it is only a matter of plugging in the ansatz into the ODE and expanding in a power series. The preceding analysis converts the equation into a finite dimensional linear algebra problem on the parameters of the $P/Q$ ansatz.
According to the suggestions from @AVK,
I built the problem in MATLAB(but I have no MATLAB to test)
I have read that example and try to adapt it to my case .
$$
\frac{y''}{(1+y'^{2})^{\frac{3}{2}}} -p - A(x-B)^2 =0
$$
given
$$
\begin{equation}
\begin{cases}
y(x_a)=YA \\
y(x_c)=YC \\
y(x_b)=YB
\end{cases}\,
\end{equation}
$$
where $ x_a \lt x_c \lt x_b$
Let $y_1=y,\quad y_2=y'$
$$
\begin{equation}
\begin{cases}
y_1'=y_2, \\
y_2'=(1+y_2^{2})^{\frac{3}{2}}\cdot (p + A(x-B)^2) =0
\end{cases}\,
\end{equation}
$$
function dydx = fcn(x,y,regin,p) % equation to solve
A = 0.001
B = 1.0e-3
dydx = zeros(2,1);
dydx = [y(2)
(1+y(2)^2)^(3/2)*(p+A*(x-B)^2)];
end
I apply the following BCs.
function res = bcfcn(ya,yb,p) % boundary conditions
YA = 0.0
YB = 5.0
YC = 3.0
res = [ya(1,1)-YA % y1(xa)=ya
ya(1,2)-YC % y1(xc)=yc Continuity
yb(1,1)-YC % y1(xc)=yc Continuity
ya(2,2)-yb(2,1) % y2(xc)=y2(xc) Continuity
yb(1,2)-YB]; % y1(xb)=yb
end
initialization
xc = 1.0;
xa = 0.0;
xb = 2.0
p = 0.0
xmesh = [0 0.25 0.5 0.75 xc xc 1.25 1.5 1.75 2];
yinit = [0.0; 0.0];
solinit = bvpinit(xmesh,yinit,p);
solve
sol = bvp5c(@fcn, @bcfcn, solinit);
fprintf('sol.parameters %7.3f.\n',...
sol.parameters);
%check solution
ypred = deval(sol,[xa,xc,xb])
fprintf("%7.3f",ypred(1,:))
I don't know if it is correct.
Because I have no MATLAB on hand, I just check it in mathworks online site.
I want to try a Fortran solver COLNEW to do similar work.
If someone has experience on this, please help me.
[ADD]
I also found a newer solver BVP_SOLVER
. Though it is easier to use, it only supports two-point BVPs and separated BCs, it is not convenient to use on my problem.
[Update]
Finally, I have solved the problem(MP-BVP with unknown parameters) using COLNEW, which is really a great solver.
Best Answer
Pedestrian answer.
In practice $h>0$ is a small number. In what follows, it will be enough that $h\in(0,1)$. Because of $$\left(e^{(\frac1h-1)x}u_1\right)'=\frac1h e^{(\frac1h-1)x}f,$$ we must have $$u_1=ce^{(\frac1h-1)x}+\int^x\frac1h e^{(\frac1h-1)(y-x)}f(y)dy.$$ There remains to determine the constant $c$, as well as the over limit point $\pm\infty$ in the integration. The key observation is that the function $$\phi_h(x):= \frac1h e^{(\frac1h-1)x}{\bf1}(x<0)$$ is integrable. A natural assumption being that the initial data $u_0(x)$ be in some $L^p$-space ($1\le p\le\infty$), we must assume as well $f\in L^p({\mathbb R})$. Then $\phi_h\star f\in L^p({\mathbb R})$, and this is the only solution belonging to $L^p({\mathbb R})$. Therefore the right choice is $$u_1=\int^x_{-\infty}\frac1h e^{(\frac1h-1)(y-x)}f(y)dy.$$ Notice that $$\|u_1\|_p\le\|\phi_h\|_1\|f\|_p=\frac1{1-h}\|f\|_p.$$ In other words $\|u^h(nh)\|_p\le(1-h)^{-n}\|u_0\|_p$. This is reminiscent to the expected $e^{-t}$ decay in the inviscid limit $h\to0+$.