I have the following nonlinear system and accompanying ODE:


My solution in MATLAB:

% Solving with ODE45
dt = 10^-1; % [s]
T = 100; % [s]
N = 2; % number of masses
X_0 = zeros(2*N,1);
t_span = [0:dt:T];
m_i = 1; % [kg]
k_i = 2; % [N/m]
c_i = 1; % [kg/s]
gamma_i = 1; % [N/m^3]
f_i = 1; % [N]
f = @(t) f_i.*sin(alpha.*t).*cos(beta.*t); % [N]
m = m_i.*ones(N,1); % [kg]
k = k_i.*ones(N+1,1); % [N/m]
c = c_i.*ones(N+1,1); % [kg/s] 
gamma = gamma_i.*ones(N+1,1); % [N/m^3]
alpha = 0.2;
beta = 0.15;
K = karray(k,gamma,N);
M = marray(m);
C = carray(c,N);
F = @(t) farray(f,t,N);
options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-3);
fun = @(t,X)odefun(X,K,C,M,F(t),N); % Defines the ODE function and it's input matrices
[t_ode,X_answer] = ode45(fun,tspan,X_0); % Calls the ODE45 function (see end of code), sets the time and axial conditions. Outputs time and displacements.

The functions karray(),carray(),marray(),farray() return array representations of the system parameters while odefun() returns an array of the EOM’s displacements and velocities. I had set the initial displacement and velocity of each mass to zero as my initial values. My confusion or question is how my kinematic constraint (i.e. zero displacement at each wall) is accounted for or enforced in the ode45() solution. Is my ode45() solution valid for the system shown or do I need to form this as a boundary value problem with the boundary conditions as x(0) = 0 and x(l) = 0 instead of an initial value problem as I have done with ode45()? Another reason I bring this up is because my ode45() solver frequently fails after t<10^-5 s when solving. Is this because there are are singularities at the walls that I am not accounting for with my ode45() solution?