Breakdown of Analytical Solution to 4th order ODE The Problem:
I have the 4th order Ordinary Differential Equation 
$$
\frac{\text{d}^4\theta}{\text{d}\eta^4}
+R(\theta-\theta_*)=0
$$
in the interval $0\le\eta\le1$, subject to the boundary conditions
$$
\eta=0: \frac{\text{d}\theta}{\text{d}\eta}=-1 ;
        \frac{\text{d}^2\theta}{\text{d}\eta^2}=0
$$
$$
\eta=1: \frac{\text{d}\theta}{\text{d}\eta}=0 ;
        \frac{\text{d}^2\theta}{\text{d}\eta^2}=0
$$
and where $\theta_*$ is to be determined such that the clamping constraint
$$\theta(\eta=0)=0$$ 
is satisfied. Skipping details, it can be shown that the solution to the differential equation is
$$
\theta=\theta_* + e^{P(\eta-1)}(A\cos P\eta+B\sin P\eta)
+e^{-P\eta}(C\cos P\eta+D\sin P\eta)
$$
where $P=\frac{R^\frac{1}{4}}{\sqrt{2}}$ and (skipping details again) the constants $A,B,C,D$ and $\theta_*$ can be determined from the boundary conditions and constraint. So far, so good.
The issue:
The solution works beautifully, until $R$ approaches $10^7$ whereupon it breaks down due to what I believe is the stiffness of the differential equation - the difference between the largest and smallest roots of the characteristic equation is of the order of $2P$ ~$R^\frac{1}{4}$. This is also apparent from the original differential equation itself, where as $R$ becomes very large $\theta \rightarrow \theta_*$ which tends to violate the Neumann boundary condition $\frac{\text{d}\theta}{\text{d}\eta}(\eta=0)=-1$. What I find very odd however, is that the breakdown in the analytical solution is manifested not at $\eta=0$, where the Neumann BC is actually satisfied very well, but by blowing up in the vicinity of $\eta=1$. This is evident in the graphic below:

My Question
Given that the analytical solution tends to break down at large $R$, how much confidence can I place in the computed values in the vicinity of $\eta=0$. The Neumann condition at $\eta=0$ certainly seems to be honoured for $R=10^7$, but I'm a bit circumspect about the correctness of the peak value in the second derivative (right plot in the graphic above). 
Any advice? Thanks in advance.
Note that in practice, I clamp the value of $\eta$ used to compute $\theta$ and its derivatives at $\eta=1.1-0.1\log_{10}R$, for $R\ge 10^6$
 A: With a proper boundary value solver there is no such boundary layer problem. Using the one from python scipy.integrate the code is
import numpy as np
from scipy.integrate import solve_bvp, odeint
import matplotlib.pyplot as plt

def odesys(t,u,R):
    th = u[0]; th_ast = u[4]; 
    return [ u[1], u[2], u[3], -R*(th-th_ast), 0*th_ast]

def boundary(u0, u1):
    return [ u0[0], u0[1]+1, u0[2], u1[1], u1[2] ]

x = np.linspace(0, 1, 3)

for R in [ 1e3, 1e4, 1e5, 1e6, 1e7 ]:
    res = solve_bvp(lambda t,u: odesys(t,u,R), boundary, x, [-x, 0*x+1, 0*x, 0*x, 0*x-1])
    x_plot = np.linspace(0, 1, 500)
    u_plot = res.sol(x_plot)
    plt.subplot(1,2,1)
    plt.plot(x_plot, u_plot[0], label='R=%.2e'%R)
    plt.subplot(1,2,2)
    plt.plot(x_plot, u_plot[2], label='R=%.2e'%R)
vlabels = [" ", "$\\theta$", "$\\theta''$"]
for k in [1,2]:
    plt.subplot(1,2,k)
    plt.legend()
    plt.xlabel("$\eta$")
    plt.ylabel(vlabels[k])
    plt.grid(); 
plt.show()

Note that $\theta_*$ was included in the state as a constant function with derivative $0$ so that the bvp solver automatically also adapts this value. This system with now 5 state dimensions allows to give all 5 boundary conditions at once.
The resulting diagram is

A: I've figured out how to avoid the issue with the breakdown of the analytical solution as follows:
From the constraint $\theta(0)=0)$ we have
$$\theta_* = -(Ae^{-P} +C)$$
and from the no slip condition $\theta''(0)=0$, we have
$$D=Be^{-P}$$
whereupon the remaining boundary conditions yield the following 3 $\times$ 3 system for the constants $A,B$ and $C$ 
$$ \begin{bmatrix}
    e^{-P} &  2e^{-P} & -1 \\
  \cos P - \sin P\ &(1+e^{-P})\cos P +(1-e^{-P})\sin P  & -e^{-P}(\cos P + \sin P)\\
-\sin P & (1-e^{-P})\cos P &e^{-P}\sin P
  \end{bmatrix}
\begin{bmatrix} A \\ B \\ C\end{bmatrix}
=
\begin{bmatrix} -P^{-1} \\ 0 \\ 0\end{bmatrix}
$$
For very large values of $R$ (large $P$) , the associated very small values of $e^{-P}$ were causing roundoff errors with the default  double precision in the VBA implementation and yielding erroneous evaluations of $B$ and $C$. Assuming a machine roundoff threshold of $10^{-16}$ (typical for double precision), this occurs at a value of $R=7.4\times 10^6$ which is just about where I noticed the analytical solution breaking down! At such low values of  $e^{-P}$, the linear system above is functionally equivalent to
$$ \begin{bmatrix}
    0 &  0 & -1 \\
  \cos P - \sin P\ &\cos P +\sin P  & 0\\
-\sin P & \cos P &0
  \end{bmatrix}
\begin{bmatrix} A \\ B \\ C\end{bmatrix}
=
\begin{bmatrix} -P^{-1} \\ 0 \\ 0\end{bmatrix}
$$
which gives $C=P^{-1}$, $A=B=0$ which yields the correct solution without any restrictions on $R$ and without recourse to the $\eta$- cutoff hack that I was employing earlier. 
So in essence, the solution to the problem that I observed is to just set $C=P^{-1}$, $A=B=0$ for $R>7.4\times 10^6$. As seen in the first graphic below, this does result in any perceptibly abrupt changes in $A$, $B$ or $C$ when this fix is implemented.

The temperature and velocity profiles now look like this

A: I found an even simpler solution - balancing the 3$\times$3 system prior to inversion for the coefficients $A,B$ and $C$.
