I have the second order boundary value problem
$$ y''=\frac{k}{2\sqrt{t}}\left[Py'-2yy'-y(P-y))\frac{1}{2t} \right], $$ with boundary conditions I will specify later in the code.
Setting $y=y_1$ and $y_2=y'$, I can consequently write the above as the following first order differential equation
$$ y_2'=y''= \frac{k}{2\sqrt{t}}\left[Py'_2-2y_1y_2-y_1(P-y_1))\frac{1}{2t} \right] $$
I am having issues implementing this on MATLAB due to the independent variable $t$. Thus far, I have written
% input parameters
P = 100000;
H0 = 18;
H1 = 3629;
t1 = 11;
k = (1/(P*sqrt(t1)))*log((H1*(P-H0))/(H0*(P-H1)));
% boundary values
ya(1) = 18;
yb(1) = 3629;
% differential equation function
f=@(x,y,t)[y(2);k/(2*sqrt(t))*(P*y(2)-2*y(1)*y(2)-y(1)*(P-y(1))/(2*t))];
bc=@(ya,yb) [ya(1)-18;yb(1)-3629]; % boundary conditions
t= linspace (0,45,100); % create mesh
solinit = bvpinit(t, [1 0]); % initial guess of the solution
sol = bvp4c(f, bc, solinit); % run solver
Each time, I get an error, however, e.g., this time that there are not enough input arguments. How should I incorporate this independent variable into the solver? If $t$ were constant, it would be easy.
If you want to get a result and are not married to the idea of just using the equation in its original form, then you can first notice that $(y(P-y))'=(P-2y)y'$, so that all the terms with $y$ on the right side are contained in a factor of the form $2tf'(y)-f(y)$. An integrating factor for that is $\frac1{2t^{3/2}}$, so that with $$ v=\frac{y(P-y)}{\sqrt{t}} $$ the differential equation gets the compact form $$ y''=\frac{k}{2}v'\implies y'=v+C=\frac{k}{2}\frac{y(P-y)}{\sqrt{t}}+C $$ This now is a Riccati equation, thus one can set $y=\frac{2\sqrt{t}}{k}\frac{u'}{u}$ to get $$ \frac{2\sqrt{t}}{k}\frac{u''}{u}-\frac{2\sqrt{t}}{k}\frac{u'^2}{u^2}+\frac{1}{k\sqrt{t}}\frac{u'}{u}=P\frac{u'}{u}-\frac{2\sqrt{t}}{k}\frac{u'^2}{u^2}+C \\~\\ 2tu''+(1-Pk\sqrt{t})u'-Ck\sqrt{t}u=0 $$ Apparently the singularity at $t=0$ still has bounded solutions that are power series in $\sqrt{t}$? Try $u(t)=f(\sqrt{t})$, or $f(s)=u(s^2)$, then $f'(s)=2su'(s^2)$, $f''(s)=4s^2u''(s^2)+2u'(s^2)=4tu''(t)+2u'(t)$ giving $$ \frac12f''(s)-\frac{Pk}{2}f'(s)-Cksf(s)=0. $$ This now is quite regular.
Let's finally translate the boundary conditions $y(0)=y_0$, $y(L)=y_L$. The Riccati transformation is singular at $t=0$, meaning that $u'(0)=\pm\infty$, which can not be implemented in a direct manner.
Thus exchange the order of the tricks, set $g(s)=y(s^2)$ in the first-order equation $$ g'(s)=2sy(s^2)=ky(P-y)+2sC=kg(P-g)+2sC $$ This is now a quite regular Riccati equation with boundary conditions $$g(0)=y_0, ~~~ g(\sqrt{L})=y_L.$$
This can now be directly implemented. In Matlab the state space for this is $(g,C)$ with $C'=0$.
One could also proceed to the Riccati transformation into a linear DE which usually serves to postpone the occurrence of singularities from the solver to the post-processing. Here one would set $g(s)-\frac{P}{2}=\frac1{k}\frac{f'(s)}{f(s)}$ to get $$ \frac{f''}{f}-\frac{f'^2}{f^2}=-\frac{f'^2}{f^2}+\frac{k^2P^2}{4}+2kCs \\~\\ f''=\left(\frac{k^2P^2}{4}+2kCs\right)f $$ with state components $(f,f',C)$, boundary conditions $$ f(0)=1/k,~~~ f'(0)=y_0-P/2,~~~ f'(\sqrt{L})-k(y_L-P/2)f(\sqrt{L})=0. $$