How to return values of differential, not the function, from odefun, ode45?

Hello this is my code: it worked fine but I need to get values of function dy2dt itself, as different column of X than values of y2 function. I tried to add y5 but not sure what write as a dy5dt equation. How could I achieve it? Other ways to get the values of dy2dt also would be fine :)

param.m1 = 3;
param.m2 = 5;
param.l0 = 5;
param.k = 25;
param.l = 10;
param.g = 10;

tspan = [0 20];
ic = [5; 1; 0; 1; 0];
[t, X] = ode45(@odeFun, tspan, ic, [], param);

function dXdt = odeFun(t, X, param)
    y1 = X(1);
    y2 = X(2);
    y3 = X(3);
    y4 = X(4);
    %y(5) = diff(y2)
    y5 = X(5);
    
    %equations k = 25
    dy1dt = y2; 
    dy2dt = (param.m2*y4^2*y1-param.k*y1+param.k*param.l0)/param.m2;
    dy3dt = y4;
    dy4dt = (-0.5*param.m1*param.g*param.l*sin(y3))/(1/3*param.m1*param.l^2+param.m2*y1^2);
    %extra equation
    dy5dt = diff(dy2dt); %error
    
    dXdt = [dy1dt; dy2dt; dy3dt; dy4dt; dy5dt];
end

Answer 1

[Jan]

You want the values of dy2dt? But you do have the equation for it already:

dy2dt = (param.m2*y4^2*y1-param.k*y1+param.k*param.l0)/param.m2;

Answer 2

[Cris LaPierre]

It would appear the values of your function dy2dt are the 2nd column of your output variable X already.

If it is diff(dy2dt) you want to return (the 5th output), the error is because, inside your function, dy2dt is a single number, so taking the diff of it results in an empty matrix.

You don't use y(d) in the solution, so if you do need to take the diff, my suggestion would be to solve your ode, then take the diff of the resulting X(:,2).

param.m1 = 3;
param.m2 = 5;
param.l0 = 5;
param.k = 25;
param.l = 10;
param.g = 10;

tspan = [0 20];
ic = [5; 1; 0; 1];
[t, X] = ode45(@odeFun, tspan, ic, [], param);

dy2dt_diff = diff(X(:,2))

dy2dt_diff = 216×1

0.0003 0.0003 0.0003 0.0003 0.0013 0.0013 0.0013 0.0013 0.0063 0.0063

●

function dXdt = odeFun(t, X, param)
    y1 = X(1);
    y2 = X(2);
    y3 = X(3);
    y4 = X(4);
    
    %equations k = 25
    dy1dt = y2; 
    dy2dt = (param.m2*y4^2*y1-param.k*y1+param.k*param.l0)/param.m2;
    dy3dt = y4;
    dy4dt = (-0.5*param.m1*param.g*param.l*sin(y3))/(1/3*param.m1*param.l^2+param.m2*y1^2);
    
    dXdt = [dy1dt; dy2dt; dy3dt; dy4dt];
end