There is no need to use inaccurate finite difference approximation. You already have an exact expression for your derivatives: the functions MyFunctionL and MyFunction. You simply need to evaluate them using the x and y values obtained from numerical integration with ode45. To do this you can either vectorize these functions or wrap them in for loops. Additionally, you should not be using global variables. That's not what they're for and they slow down your code. You should pass additional parameter via an anonymous function. Here's how you might re-write your two integration functions to vectorize them and not use globals:
function dy = MyFunctionL(x,y,Fg,Lg,EI)
dy(1,:) = y(2,:);
dy(2,:) = (Fg/EI)*(Lg-x);
function dy = MyFunction(x,y,Fg,Lg,EI)
dy(1,:) = y(2,:);
dy(2,:) = (Fg/EI)*(Lg-x).*(1+y(2,:).^2).^1.5;
You can then call the sub-functions (they could be in separate M-files too), passing in the additional parameters, and obtain the derivatives like this:
[x,y] = ode45(@(x,y)MyFunctionL(x,y,Fg,Lg,EI),xspan,x0);
xdot = MyFunctionL(x.',y.',Fg,Lg,EI).';
...
[x,z] = ode45(@(x,y)MyFunction(x,y,Fg,Lg,EI),xspan,x0);
zdot = MyFunctionL(x.',y.',Fg,Lg,EI).';
Note the transposes (.') which are present because ode45 needs a column vector but outputs each time step along the columns (for performance reasons). You could equally write separate functions that use your expressions for the derivates.
Analytic solution :
$$y(x)=-\int \frac{-ax+b+k}{-x^2+x+m} = \frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m | +C$$
With condition $y(0)=-y(1)$ :
$y(0)=\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m | +C$
$C=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|$
$$y(x)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$
As a consequence :
$y(1)=y(0)-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m| + {\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|m|} $
$y(1)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
The second condition $y(1)=-y(0)$ implies :
$-y(0)=y(0)+2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
$y(0)=-\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{1}{\sqrt{4m+1}} \right)$
$$y(x)=-2\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{-1}{\sqrt{4m+1}} \right) + \frac{a}{2}\ln|m|+ +{\frac{a-2b-2k}{\sqrt{4m+1} }\tanh^{-1}\left( \frac{2x-1}{\sqrt{4m+1}} \right) - \frac{a}{2}\ln|-x^2+x+m |} $$
Best Answer
ode45 is for initial value problems: all conditions occur at the same value of the independent variable $x$. It looks to me like what you want is a boundary value problem, where some conditions are at one endpoint of the interval and others at the other endpoint. That's what bvp4c is for.