4.6 Solve homogeneous 2nd order ODE, constant coefficients (BVP)

Problem: Solve

\[ y^{\prime \prime }\left ( t\right ) +t\ y\left ( t\right ) =0 \]

with the boundary conditions

\[ y\left ( 0\right ) =3,y\left ( 20\right ) =-1 \]

For solving with Matlab, the ODE is first converted to state space as follows

Given \(y^{\prime \prime }\left ( t\right ) +t\ y\left ( t\right ) =0\), let

\begin{align*} x_{1} & =y\\ x_{2} & =y^{\prime }\\ & =x_{1}^{\prime }\end{align*}

Therefore

\[ x_{1}^{\prime }=x_{2}\]

And

\begin{align*} x_{2}^{\prime } & =y^{\prime \prime }\\ & =-t\ y\\ & =-t\ x_{1}\end{align*}

This results in

\[\begin {bmatrix} x_{1}^{\prime }\\ x_{2}^{\prime }\end {bmatrix} =\begin {bmatrix} 0 & 1\\ -t & 0 \end {bmatrix}\begin {bmatrix} x_{1}\\ x_{2}\end {bmatrix} \]

Now bvp4c() can be used.

Mathematica

Clear[y,t]; 
eq  = y''[t]+t y[t]==0; 
ic  = {y[0]==3,y[20]==-1}; 
sol = First@DSolve[{eq,ic},y[t],t]; 
 
y   = y[t]/.sol
 

\[ \frac {-\sqrt {3} \text {Ai}\left (\sqrt [3]{-1} t\right )+\text {Bi}\left (\sqrt [3]{-1} t\right )-3\ 3^{2/3} \Gamma \left (\frac {2}{3}\right ) \text {Bi}\left (20 \sqrt [3]{-1}\right ) \text {Ai}\left (\sqrt [3]{-1} t\right )+3\ 3^{2/3} \Gamma \left (\frac {2}{3}\right ) \text {Ai}\left (20 \sqrt [3]{-1}\right ) \text {Bi}\left (\sqrt [3]{-1} t\right )}{\sqrt {3} \text {Ai}\left (20 \sqrt [3]{-1}\right )-\text {Bi}\left (20 \sqrt [3]{-1}\right )} \]

Plot[Re[y],{t,0,20}, 
     FrameLabel->{{"y(t)",None}, 
      {"t","Solution"}}, 
     Frame->True, 
     GridLines->Automatic, 
     GridLinesStyle->Automatic, 
     RotateLabel->False, 
     ImageSize->300, 
     AspectRatio->1, 
     PlotRange->All, 
     PlotStyle->{Thick,Red}, 
     Exclusions->None]
 

Matlab

function e56 
 
t0 = 0; %initial time 
tf = 20; %final time 
 
solinit = bvpinit(linspace(t0,tf,5),[3 -1]); 
 
sol = bvp4c(@odefun,@bcfun,solinit); 
 
plot(sol.x(1,:),sol.y(1,:),'r') 
title('solution'); 
xlabel('time'); 
ylabel('y(t)'); 
grid; 
set(gcf,'Position',[10,10,320,320]); 
 
function dydt=odefun(t,y) 
dydt=[y(2); -t.*y(1)]; 
 
function res=bcfun(yleft,yright) 
res=[yleft(1)-3 
    yright(1)+1];