How to solve following BVP ODE?

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".

I am looking forward to hearing your advices.

Cheers,

Dariush

Answer 1

[Alan Stevens]

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]

Thanks Alan for the answer! It is exactly what I was trying to do.

cheers,

Dariush