How to plot not only the result but also the derivatives of an ode using the ode45 function in Matlab? I have already successfully run a code for the simulation of the deflection of beams under different loadings. I used the Matlab program and the ode45 solver for Initial value Problems and bvp4c solver for Boundary conditions. So I generated graphics such as the one below:

I want now to plot not only the deflection y(x) but also dy/dx and d^2y/dx^2 on the same graphic and I am having difficulties to actually do it. Hope someone can give me a hint. All I want is to plot the function I inserted on the MyFunctionL and MyFunctionNL and it´s derivative together with the result. One of the codes I generated is below:
function def1= def1(F,L)
%--------------------------------------------------------------
% Deflection 
% F [N] und L [mm]
% EI N*mm^2
% y[mm]
%--------------------------------------------------------------
global Fg Lg EI;
Fg = F;
Lg = L;
Em=200*10^9

%This part below is not relevant for the question
 for i=1:3
      if i==1
          b=0.055
          h=0.1
          Im=b*h^3/12
      end
      if i==2
          Im=2.5175*10^(-6)
      end
      if i==3
          re=0.065
          ri=0.060
          Im=pi/4*((re^4)-(ri^4))
      end

%As Im is in m^4 we are converting EI to N*mm^2 by multiplying it by (1000^2).    
EI=(Em*Im)*(1000^2)

$Now this part below is.
%Start point        
x0 = [0 0];
%Längenintervall
xspan = [0 L];
%Call Solver -> Linear
[x y] = ode45(@MyFunctionL,xspan, x0);
%Plot result
figure
plot(x,y(:,1));
if i==1  title('y(x) in a Retangular Profile') 
end
if i==2  title('y(x) in a Iprofile(IPE 100)')  
end
if i==3  title('y(x) in a Round hollow section') 
end
set(gca,'ydir','reverse')
xlabel('[mm]');
ylabel('[mm]');
hold on;
%Call Solver -> NonLinear
[x z] = ode45(@MyFunction,xspan, x0);
%Plot result
plot(x,z(:,1),'-r');
set(gca,'ydir','reverse')
end   
return

%---------------------------------------------------------------
function dy = MyFunctionL(x,y)
global Fg Lg EI;
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = (Fg/EI)*(Lg - x);
return

function dy = MyFunction(x,y)
global Fg Lg EI;
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = (Fg/EI)*(Lg - x)*(1+y(2)^2)^1.5;
return

 A: In the string

[x y] = ode45(@MyFunctionL,xspan, x0);

the function ode45 returns not only $y(1)$ (this is your $y(x)$), but also  $y(2)$ (this is you $\frac{dy}{dx}$). Therefore, to plot both of them on the one figure, you just need to write

plot(x,y);

instead of

plot(x,y(:,1))).

But note that the value of your derivative is very small in comparison with $y(x)$. (It will be slightly difficult to see - at the top of your plot).
However, the graph of $\frac{d^2y}{dx^2}$ is a bit more tricky. You need to use some difference approximation, for example like this:

YY = diff(Y(:,2)) ./ diff(X);
YY(end+1) = 0; %it is the technical string
YY = [Y YY];
plot(X, YY);

Of course, this way maybe not optimal. But at least something.
A: 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.
