Solve second order differential equation with multidimensional time-dependent matrix as coefficients?

I am trying to solve a multidimensional mass-damper system with time-dependent mass coefficient and damping coefficient. The differential equation describing the model is as following:

I have tried solving it with the usual ode45 solver. As parameter, I call in the funtion ParaFunct, where in its function it calls the function callMat. callMat is the function that calls a certain coefficient value (this would be , ), and , dependent of the time t. However, the data matrices from which I call in my time-dependent parameters is just a matrix of discrete values (MatrixStruct). It is only possible to fit it with pchip or spline. In conclusion, I think the current parameter calling funtion is a bad idea, because of the discrete values.

My question is: Say I fitted my time dependent parameters with either pchip or spline. Is is then a good idea to use matrix of function handles while using ode45 solver? if no, what are the better alternatives to my problem?

%%ParaFunct
function xp = paraFunc(t, x, MatrixStruct,ForceMat, t_ind, mass)
xp=zeros(2,1);
[A, Force]= callMat(MatrixStruct, ForceMat, t, t_ind);
% xp(1)=x(2);
%xp(2)=-D_mat*x(2)-K_mat*x(1)+Force;
B= zeros(size(A,1),1);
B(size(A,1)/2+1,1)=1/mass;
xp=A*x+B.*Force;

%%callMat
function [A,Force]=callMat(MatrixStruct, ForceMat, t, t_ind)
[o,p]=min(abs(t_ind-t));
p
A=MatrixStruct(p).mat;
Force=ForceMat(p);
end

%%odesolver
Force= 100;
t_ind= 0:0.1:4.6;
ForceMat=[Force*ones(1,40), zeros(1,40)];
    tspan=[min(t_ind) max(t_ind)];
    
    options = odeset('AbsTol',1e-9,'RelTol',1e-9);
    [t,y]=ode45(@(t,y) paraFunc(t, y, MatrixStruct,ForceMat, t_ind, mass), tspan, x, options);