Modelling Equations of Motion using ODE45

So I want to model equation 8 (Equation of motions) and then run ODE45 to solve it. I thought this might involve converting 2nd order to 1st order by doing the following:

and code it like this:

but i might need some more help

%M*x‘‘(t)+D*x‘(t)+K*x(t)= F(t)
M=1;
omegax=0.9;
omegaz=1.1;
omega=0.5;
alpha=0.1;
gammax=0.02;
gammaz=0.04;
mu=0.41;

function main
x0=[0 1];      % first component is x(0), second x'(0)
tspan=[0 20];  
[T,Y] = ode45(@myfun,tspan,x0);
plot(T,Y(:,1),T,Y(:,2));
function dy = myfun(t,y)
F=??????????
K=[omegax^2, -omega^2; -omega^2, omegaz^2];
D=[gammax, 0; 0 ,gammaz]
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=(F-K*y(1)-D*y(2))/M; %M*x‘‘(t)+D*x‘(t)+K*x(t)= F(t) rearranged for x''

[Torsten]

You need 4 equations to define your system, not 2: x,x',z,z' are to be solved for.

[William Ta]

How would the 4 equations look like?

[Bjorn Gustavsson]

It is as simple as:

function drdtdvdt = eq_of_motion(t,rv,more,input-arguments,and,such)
drdt = rv(3:4); % 
F = your calculation of the force
dvdt = F/m;
drdtdvdt = [drdt(:);dvdt(:)];

Position, velocity and acceleration will be a 2x1 array in your case, therefore you need to keep track of both components for each, and you're almost there.

HTH

[Huda]

How I can write F=?????? in MATLAB

[Sam Chak]

Hi @Huda

If you have a similar problem that requires determing the "F" or the Force thing, can you post a new question?

The F is not given in this problem, but there might other ways to mathematically describe F or to numerically produce the data points "F" subject to certain cost functions and constraints, without destabilizing the system.