# How to solve following BVP ODE?

3 views (last 30 days)
Dariush Javani on 13 May 2021
Commented: Dariush Javani on 13 May 2021
Hi mates,
I am trying to solve the following ODE in MATLAB, but I do not attain a reasonable answer. It would be highly appreciated if you let me know of the bugs in my code below the ODE.
The ODE is: (d/dx)(y^3 dy/dx)+(2/3)(x dy/dx)-(1/3)y=0;
The BC's are: (y^3)(dy/dx)=-1, for x=0; and y(inf)=0;
What I have done is to use transforms as below
Y(1)=y;
Y(2)=dy/dx;
Then, dY(1)/dx=Y(2) and dY(2)/dx=((1/3)(1/Y(1)^2))-((2/3)(xY(2)/Y(1)^3)) (to have this equation I just ignored the high order small value of (dy/dx)^2 which is an assumption I took; If there is better solution I would be happy to replace).
According to above descriptions, I used following code in MATLAB:
function liftoff()
xmesh=linspace(0,4,10);
solinit = bvpinit(xmesh,[1 0]);
sol = bvp4c(@lode,@bcs,solinit);
plot(sol.x,sol.y(1,:),'b-x');
function dYdx = lode(x,Y)
dYdx(1) = Y(2);
dYdx(2) = (1/3)*((1/(Y(1)^2))-(2*x*Y(2)/(Y(1)^3)));
function res = bcs(ya,yb)
res = [ ((ya(1)^3)*ya(2))+1;
yb(1)];
But, I encountered with error of "a singular Jacobian".
Cheers,
Dariush

Alan Stevens on 13 May 2021
Weird ODE and initial boundary condition! You can solve for small values of x using the following (where y(x=0) was chosen by trying a few values). Larger values of y0 lead to an upturn in y after the dive!
y0 = 1.311435;
v0 = -1./y0.^3;
Y0 = [y0, v0];
xend = 4;
xspan = [0 xend];
[x, Y] = ode15s(@fn, xspan, Y0);
y = Y(:,1);
v = Y(:,2);
plot(x,y),grid
function dYdx = fn(x,Y)
y = Y(1);
v = Y(2);
if y<=0 % just to prevent weird behaviour !!
dYdx = [0; 0];
else
dYdx = [v;
(y/3 - 2*x.*v/3 - 3*y.^2.*v.^2 )./y.^3];
end
end
Dariush Javani on 13 May 2021
Thanks Alan for the answer! It is exactly what I was trying to do.
cheers,
Dariush

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!