How to plot not only the result but also the derivatives using the ode45 function?

Hello there, I have already successfully run a code for the simulation of the deflection of beams under different loadings. I used the ode45 for Initial value Problems and bvp4c 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

Answer 1

[Star Strider]

The ODE solvers return (in your situation) [y(x) y’(x)] and [z(x) z’(x)], so the first derivative is in y(:,2) and z(:,2). I didn’t run your code, but I ran a similar version of these functions with an ODE I have been working with, and since this strategy worked with it, I adapted yours to:

MyDrvsL = @(x,y) [y(:,2)';  (Fg/EI).*(Lg - x')];
MyDrvs  = @(x,y) [y(:,2)';  (Fg/EI).*(Lg - x').*(1+y(:,2)'.^2).^1.5];

Given (Nx2) vectors for y and z, these should return [y’ y”]’ as a (2xN) matrix. (Note the transposition.)

[Daniel Clemente]

Thanks!

[Star Strider]

My pleasure!