multiplying a function handle by a constant

For the code below I am trying to find the polynomial equation that represents the system. There are 4-2nd ODE equations I have made them into 4-1st ODE the first order differentials become "y1,y2,y3and y4" and the 2nd order differentials become a first order. I then put all 4 into a matlabFunction, how do I multiple the function handle by the constant"h" with out getting the error "Operator ‘*’ is not supported for operands of type ‘function_handle’."?
The rest of the code works if working with a single differential equation fx(x,y,t) with only "x,y,t" but in my equations I have "Xf,Xr,Xb,theta" and "Vf,Vr,Vb,omega" that I have chosen to represent as "x1,x2,x3,x4" and "y1,y2,y3,y4" respectively. The next question is will I run into problems here as well?
I am not that expirenced working with matlabFunction command to know how to get this to work. The project requires me to use 2 different numerical methods to find the polynomial equation that best fits so I cannot use "Polyfit" to get the polynomial that best fits. Any suggestions that can help me to get this to work would be appreciated.
clear,close,clc

%______________________________________________________________SOLUITION_2_Heun’s_Method_(for second order differential equations)_&__Least-Square_Nethod____________________________________________________________%

%4 Equations representing the system working with

% MfXf"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb’-(L1*theta’)-Xf’)-(Kf*Xf);
% MrXr"=Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb’+(L2*theta’)]-Xr’)-(Kr*Xr) ;
% MbXb"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb’-(L1*theta’)]-Xf’)+Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb’+(L2*theta’)]-Xr’)+fa(t);
% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf’-(L1*theta’))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr’+(L2*theta’))*L2]+[fa(t)*L3]};
clc
clear
%————————————————-SYSTEM_PARAMETERS———————————————————————————————————————————————————-
Ic=1356; %kg-m^2
Mb=730; %kg
Mf=59; %kg
Mr=45; %kg
Kf=23000; %N/m
Ksf=18750; %N/m
Kr=16182; %N/m
Ksr=12574; %N/m
Bsf=100; %N*s/m
Bsr=100; %N*s/m

L1=1.45; %m
L2=1.39; %m
L3=0.67; %m
t=[0:20]; % time from 0 to 20 seconds
n=4; %order of the polynomial

%Initial Conditions-system at rest therefore x(0)=0 dXf/dt(0)=0 dXr/dt(0)=0 dXb/dt(0)=0 dtheta/dt(0)=0 ;
%time from 0 to 20 h=dx=5;
x0=0; %x at initial condition
y0=0; %y at initial condition
t0=0; %t at the start
dx=5; %delta(x) or h
h=dx;
tm=20; %what value of (x) you are ending at

Xf = sym(‘x1’); %x1=Xf
Xr = sym(‘x2’); %x2=Xr
Xb = sym(‘x3’); %x3=Xb
theta = sym(‘x4’); %x4=theta

Vf = sym(‘y1′); %y1=Xf’ = Vf = dXf/dt = Xf_1
Vr = sym(‘y2′); %y2=Xr’ = Vr = dXr/dt = Xr_1
Vb = sym(‘y3′); %y3=Xb’ = Vb = dXb/dt = Xb_1
omega = sym(‘y4’); %y4=theta’= omega = dtheta/dt = theta_1
t = sym(‘t’);

% MfXf"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb’-(L1*theta’)-Xf’)-(Kf*Xf);
% Vf’=(1/Mf)(Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb’-(L1*theta’)-Xf’)-(Kf*Xf));

% MrXr"=Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb’+(L2*theta’)]-Xr’)-(Kr*Xr) ;
% Vr’=(1/Mr)(Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb’+(L2*theta’)]-Xr’)-(Kr*Xr)) ;

% MbXb"=Ksf([Xb-(L1*theta)]-Xf)+Bsf([Xb’-(L1*theta’)]-Xf’)+Ksr([Xb+(L2*theta)]-Xr)+Bsr([Xb’+(L2*theta’)]-Xr’)+fa(t);
% Vb’=(1/Mb)(-Ksf([Xb-(L1*theta)]-Xf)-Bsf([Xb’-(L1*theta’)]-Xf’)-Ksr([Xb+(L2*theta)]-Xr)-Bsr([Xb’+(L2*theta’)]-Xr’)+fa(t));

% Ic*theta"={-[Ksf(Xf-(L1*theta))*L1]-[Bsf(Xf’-(L1*theta’))*L1]+[Ksr(Xr+(L2*theta))*L2]+[Bsr(Xr’+(L2*theta’))*L2]+[fa(t)*L3]};
% omega’=(1/Ic){[-Ksf*Xf + Ksf((L1)^2)*theta)] + [-Bsf*Vf*L1 + Bsf*Vf*L1 + Bsf((L1)^2)*omega)] + [Ksr*Xr*L2 + Ksr((L2)^2)*theta)] + [Bsr*Vr*L2 + Bsr((L2)^2)*omega)*L2] + [fa(t)*L3]};

Xf_2=(Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*omega)-Vf)-(Kf*Xf)))/Mf;
Xr_2=(Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)-(Kr*Xr))/Mr;
Xb_2=(Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Vb-(L1*theta’))-Vf)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Vb+(L2*omega))-Vr)+(10*exp(-(5*t))))/Mb;
theta_2=((-(Ksf*(Xf-(L1^2*theta))*L1)-(Bsf*(Vf-(L1*omega))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Vr+(L2*omega))*L2)+((10*exp(-(5*t)))*L3)))/Ic;

Eqns=[Xf_2; Xr_2; Xb_2; theta_2];
F1=matlabFunction(Eqns,’Vars’,{‘x1′,’x2′,’x3′,’x4′,’y1′,’y2′,’y3′,’y4′,’t’})

%==INPUT SECTION for Euler’s and Heun’s==%
fx=@(x,y,t)y;
fy=@(x,y,t)F1;

%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn(1) = x0;
yn(1) = y0;

for i=1:length(tn)
%==EULER’S METHOD
xn(i+1)=xn(i)+fx(xn(i),yn(i),tn(i))*h;
yn(i+1)=yn(i)+fy(xn(i),yn(i),tn(i))*h;

%==NEXT 3 LINES ARE FOR HEUN’S METHOD

tn(i+1)=tn(i)+h;

xn(i+1)=xn(i)+0.5*(fx(xn(i),yn(i),tn(i))+fx(xn(i+1),yn(i+1),tn(i+1)))*h;
yn(i+1)=yn(i)+0.5*(fy(xn(i),yn(i),tn(i))+fy(xn(i+1),yn(i+1),tn(i+1)))*h;
fprintf(‘t=%0.2ft x=%0.3ft y=%0.3fn’,tn(i),xn(i),yn(i))
end

%%%LEAST SQUARE METHOD-FINDS POLYNOMIAL FOR GIVEN DATA SET%%%%%
%INPUT SECTION for Least-Square
X=xn;
Y=yn;

%%__CALCULATIONS SECTION__%%
k=length(X); %NUMBER OF AVAILABLE DATA POINTS
m=n+1; %SIZE OF THE COEFFICENT MATRIX
A=zeros(m,m); %COEFFICENT MATRIX
for j=1:m
for i=1:m
A(j,i)=sum(X.^(i+j-2));
end
end
B=zeros(m,1); %FORCING FUNCTION VECTOR
for i=1:m;
B(i)=sum(Y.*X.^(i-1));
end
a1=AB %COEFFICIENTS FOR THE POLYNOMINAL–> y=a0+a1*x+a2*x^2….an*x^n CAN BE REPLACED BY GAUSSIAN ELIMINATION

%%%%%=========GAUSSIAN ELIMINATION TO FIND "a"========%%%%%%
%%%INPUT SECTION
%CALCULATION SECTION
AB=[A B]; %Augumentent matrix
R=size(AB,1); %# OF ROWS IN AB
C=size(AB,2); %# OF COLUMNS IN AB

%%%%FOWARD ELIMINATION SECTION
for J=1:R-1
[M,I]=max(abs(AB(J:R,J))); %M=MAXIMUM VALUE, I=LOCATION OF THE MAXIMUM VALUE IN THE 1ST ROW
temp=AB(J,:);
AB(J,:)=AB(I+(J-1),:);
AB(I+(J-1),:)=temp;

for i=(J+1):R;
if AB(i,J)~=0;
AB(i,:)=AB(i,:)-(AB(i,J)/AB(J,J))*AB(J,:);
end
end
end

%%%%BACKWARDS SUBSTITUTION
a(R)=AB(R,C)/AB(R,R);
for i=R-1:-1:1
a(i)=(AB(i,C)-AB(i,i+1:R)*a(i+1:R)’)/AB(i,i);
end
disp(a)
syms X
P=0;
for i=1:m;
TT=a(i)*X^(i-1); %T=INDIVIDUAL POLYNOMIAL TERMS
P=P+TT;
end
display(P)
%========END OF GAUSSIAN ELIMINATION=======%%%%%%%%For the code below I am trying to find the polynomial equation that represents the system. There are 4-2nd ODE equations I have made them into 4-1st ODE the first order differentials become "y1,y2,y3and y4" and the 2nd order differentials become a first order. I then put all 4 into a matlabFunction, how do I multiple the function handle by the constant"h" with out getting the error "Operator ‘*’ is not supported for operands of type ‘function_handle’."?
The rest of the code works if working with a single differential equation fx(x,y,t) with only "x,y,t" but in my equations I have "Xf,Xr,Xb,theta" and "Vf,Vr,Vb,omega" that I have chosen to represent as "x1,x2,x3,x4" and "y1,y2,y3,y4" respectively. The next question is will I run into problems here as well?
I am not that expirenced working with matlabFunction command to know how to get this to work. The project requires me to use 2 different numerical methods to find the polynomial equation that best fits so I cannot use "Polyfit" to get the polynomial that best fits. Any suggestions that can help me to get this to work would be appreciated.
clear,close,clc