How to plot not only the solution but also all derivatives using ode45 in Matlab?

Question

I have already successfully run a code for the simulation of the deflection of beams under different loading. I used Matlab and the ode45 solver for Initial value Problems.

I want now to plot not only the deflection y(x) but also dy/dx and d^2y/dx^2 on the same figure 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

Answer 1

You can use Matlab's diff() function to calculate numerical derivatives. The following assumes that the biggest value in 'x' is the last value, and that the vector is equally-spaced.

dx=x(end)/length(x);         % Calculate dx

dy=diff(y)/dx;            % Calculate numerical derivative

dx2=x(end-1)/length(dy);  % Calculate dx2

d2y=diff(dy)/dx2;      % Calculate numerical second derivative

x=x(1:end-2);                % Shorten vectors for plotting
y=y(1:end-2);
dy=dy(1:end-1);

figure;plot(x,y,x,dy,x,d2y);

If ode45() is not giving you equally-spaced time points, you should make your tspan vector like this, say if you want 1,000 points:

tspan = linspace(0,L,1000);

Also you don't need the semicolon after declaring global variables. I think that might even mess something up.

Hope that is helpful!