how to solve a higher order boundary conditioned differential equation with matlab function

I'm trying to solve a boundary value condition problem with matlab bvp4c function. But I get warning and it doesn't give correct solution.

my equation is d4y/dx4=C/(ax^2 -y)^2 ,C and a are constant , my warning content is following below

Warning: Unable to meet the tolerance without using more than 2500 mesh points.

The last mesh of 2152 points and the solution are available in the output argument.

The maximum residual is 9.25647e+11, while requested accuracy is 0.001.

my code:

clc;clear
    solinit = bvpinit([0,200],[1,1,1,1]);
    sol = bvp4c(@deriv,@bcs,solinit);
    plot(sol.x,sol.y)
    function dydx = deriv(x,y)
    
    c=1;a=1;
    dydx = [y(2)
            y(3)
            y(4)
           (c/(a*x.^2 -y(1)).^2)];
    end
    
    function res = bcs(ya,yb)
    res = [ ya(1) 
            ya(2)
            yb(3)
            yb(4)];
    end
    

I don't know what is its reason?!

thanks in advance for your helping

Answer 1

[Star Strider]

This now runs.

One problem was that ‘deriv’ is actually a reserved name, so I changed it. Also the ‘guess’ function was incorrectly coded, and ‘xmesh’ needs to be a vector of more than two elements.

Try this —

xmesh = linspace(0, 200, 5000);
guess = @(x) [(x+1); -(x+1); (2-x); (2+x)];
solinit = bvpinit(xmesh,guess);
sol = bvp4c(@bvpderiv, @bcs,solinit);

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 5000 points and the solution are available in the output argument.
The maximum residual is 343.374, while requested accuracy is 0.001.

figure
plot(sol.x,sol.y)
grid
legend(string(Subs), 'Location','best')

function dydx = bvpderiv(x,y)
c=1;a=1;
dydx = [y(2)
    y(3)
    y(4)
    (c/(a*x.^2 -y(1)).^2)];
end

function res = bcs(ya,yb)
res = [ ya(1) 
    ya(2)
    yb(3)
    yb(4)];
end

If the eventual desire is to have ‘a’ and ‘C’ as arguments to ‘bvpderiv’, rather than hard-coding them, code it as —

function dydx = bvpderiv(x,y,a,c)

dydx = [y(2)

y(3)

y(4)

(c/(a*x.^2 -y(1)).^2)];

end

and then call it as —

@(x,y)bvpderiv(x,y,a,c)

with the parameters previously defined in the workspace.

.