See the following code. The details can be quite involved, so here I just give an outline.
1. First you need to convert your 3rd order ODE to a system of 3 first order ODEs. See here for example: http://mathonline.wikidot.com/converting-nth-order-odes-to-systems-of-n-first-order-odes.
2. Then you need to write that system of ODEs into a function. See this example: https://www.mathworks.com/help/releases/R2020a/matlab/ref/bvp4c.html#mw_9083ab4d-6cde-4491-95b6-cd9b6313a459
3. Then you need to create the function to specify the boundary condition. See link in point 2 for detail.
4. Then create the initial guess. Also see link in point 2.
xmesh = linspace(0, 1, 100);
solInit = bvpinit(xmesh, [0; 0; 0]);
sol = bvp4c(@odeFun, @bcFun, solInit);
plot(sol.x, sol.y)
legend({'$y$', '$\dot{y}$', '$\ddot{y}$'}, ...
'FontSize', 16, ...
'Location', 'best', ...
'Interpreter', 'latex')
function dydx = odeFun(x, y) % ODE function
dydx = zeros(3,1);
dydx(1) = y(2);
dydx(2) = y(3);
dydx(3) = -2*exp(-3*y(1)) + 4*(1+x).^-3;
end
function res = bcFun(ya, yb) % boundary Condition fcn
res = [ya(1)-0;
ya(2)-0;
yb(1)-log(2)];
end
