Assuming your parameters are constants, create the dynamics function as follows:
function dY = DinFun(~,Y,parameters)
% constants
P = parameters.P;
m1 = parameters.m1;
m2 = parameters.m2;
R1 = parameters.R1;
R2 = parameters.R2;
K1 = parameters.K1;
K2 = parameters.K2;
% states
x1 = Y(1);
x2 = Y(2);
v1 = Y(3);
v2 = Y(4);
% state derivatives
dY(1) = v1;
dY(2) = v2;
dY(3) = (1/m1) * (P - R1*(v1 - v2) - K1*(x1-x2));
dY(4) = (1/m2) * (-R2*v2 - K2*x2 + R1*(v1-v2) + K1*(x1-x2));
dY = dY';
end
than you need a script like the following:
clc; close all; clear
%% defining the parameters involved
parameters.P = 1;
parameters.m1 = 10;
parameters.m2 = 1;
parameters.R1 = 0.1;
parameters.R2 = 0.3;
parameters.K1 = 100;
parameters.K2 = 80;
%% Initial state
% Assuming you want to start from 0,0 with null velocities
x1_0 = 0;
x2_0 = 0;
v1_0 = 0;
v2_0 = 0;
Y_0 = [x1_0; x2_0; v1_0; v2_0];
%% ode settings
% assuming your time starts from 0
t0 = 0;
tf = 10;
%% ODE
[T,Y] = ode45(@DinFun, [t0 tf], Y_0, [], parameters);
%% What ever else you need
You just need to change the value of the constants and the final time
