Here is one approach, see https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html for further details on passing extra parameters.
This example is just intended to give you an idea of how to approach this. You will have to modify it based on your needs, and actual values
% Make up some values just to demonstrate program
m = 1;
I = 2;
a = 8;
b = 5;
u = 4;
nPts = 100;
tFinal = 23;
ty = linspace(0,tFinal,nPts);
Fyf = randn(nPts,1); % time varying force Fyf
Fyr = randn(nPts,1); % time varying force Fyr
v0 = 0; % initial velocity
r0 = 0; % initial position
[t,v,r] = mySolve(m,I,a,b,u,ty,Fyf,Fyr,v0,r0);
plot(t,v,t,r);
legend('v','r')
ylabel('position and velocity')
xlabel('time')
% Everything below this line could be put into a separate file
function [t,v,r] = mySolve(m,I,a,b,u,ty,Fyf,Fyr,v0,r0)
%motionSolve - solves my system of ode's
% m,I,a,b,u - System parameters (assumed to be scalar constants)
% ty - vector of time values where forces are defined
% Fyf,Fyr - vectors of force values occuring at times corresponding to
% elements of ty
% Find start and end time
tspan = [min(ty),max(ty)];
% Set initial condition
x0 = [v0,r0];
% Solve ode
[t,x] = ode45(@odefun,tspan,x0);
% Return desired values
v = x(:,1);
r = x(:,2);
% Use nested function so that external inputs, Fyf and Fyr and parameters
% are in scope when computing derivative values
function dxdt = odefun(t,x)
% Calculates derivatives with respect to time
% Interpolate to find force values at current time
ff = interp1(ty,Fyf,t);
fr = interp1(ty,Fyr,t);
% Assign local variable names for readability
v = x(1);
r = x(2);
% Calculate derivatives
dxdt(1,1) = 1/m*(ff + fr ) - u*r;
dxdt(2,1) = 1/I*(a*ff - b*fr);
end
end
