Solving a system of second order differential equations (coupled mass-spring system)

Hello, I am trying to solve the system defined by the following equations EqX1 and EqX2:

syms X1(t) X2(t) K f m1 m2 Y t
EqX1 = diff(X1,2) == K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f + K*(X2-X1)/m2;
[ODE,Vars] = odeToVectorField(EqX1, EqX2);
Sys = matlabFunction(ODE, 'Vars', {t, Y, K, f, m1, m2});

With arbitrary constants m1, m2, f, and K.

m1 = 1;                 % Mass 1
m2 = 2;                 % Mass 2
K = 100;                % Spring constant
f = 5;                  % Forcing term on mass 2

All of my initial conditions are zero, except for the velocity of mass 2 (v2_i).

v2_i = 3;
ic = [0 v2_i 0 0];
tspan = [0, 1];
[T,Y] = ode45(@(t,Y) Sys, tspan, ic);

I feel like this should work, but when I try to run it, I get the following errors:

Error using odearguments (line 95)
@(T,Y)SYS returns a vector of length 1, but the length of initial conditions vector is 4. The vector returned by @(T,Y)SYS and the initial conditions vector must have the same
number of elements.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in CoupledMassSim (line 29)
[T,Y] = ode45(@(t,Y) Sys, tspan, ic);

Having looked at similar problems on this forum, I'm not sure why @(T,Y)SYS is returning a vector of length 1. Here is the Sys function for reference:

Sys =

  function_handle with value:

    @(t,Y,K,f,m1,m2)[Y(2);f+(K.*(Y(1)-Y(3)))./m2;Y(4);-(K.*(Y(1)-Y(3)))./m1]

Answer 1

[Stephan]

syms X1(t) X2(t) K f m1 m2 Y t
EqX1 = diff(X1,2) == K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f + K*(X2-X1)/m2;
[ODE,Vars] = odeToVectorField(EqX1, EqX2)
Sys = matlabFunction(ODE, 'Vars', {t, Y, K, f, m1, m2})
m1 = 1;                 % Mass 1
m2 = 2;                 % Mass 2
K = 100;                % Spring constant
f = 5;                  % Forcing term on mass 2
v2_i = 3;
ic = [0 v2_i 0 0]';
tspan = [0, 1];
[T,Y] = ode45(@(t,Y)Sys(t,Y,K,f,m1,m2), tspan, ic);
plot(T,Y)

[Benjamin Tonita]

Thank you!

Also turns out I had my signs turned around. Equations should be:

EqX1 = diff(X1,2) == -1*K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f - K*(X2-X1)/m2;