obtaining single polynomial equation using Heun's method from 4 second order differential equations

I get this error when trying to solve these 4 second order differential equations using a program that solves them using Heun’s method and then Gaussian elimination to get a polynomial equation the represents the system. Previsouly posted for help, but the soluitions have symbolics and I need a singular soluition that represents all 4 equations. Does any one have a suggestion on how I can slove this? Is using the "function" with all 4 equations the best way to do this, or is there a nicer method for doing this?
clear,close,clc

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

%4 Equations representing the system working with

%————————————————-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
t0=0; %time at the start
tm=20; %what value of (t time) you are ending at
t=t0:tm; % time from time start to time finish 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
dx=5; %delta(x) or h
h=dx;

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

%==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 y=%0.3ft z=%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=======%%%%%%%%

%STANDARD DEVIATION
Y_bar=mean(Y); %ADVERAGE OF y
St=sum((Y-Y_bar).^2);
SD=sqrt(St/(k-1)); %STANDARD DEVIATION

%STANDARD ERROR
for i=1:m;
T(:,i)=a(i)*X.^(i-1); %T=INDIVIDUAL POLYNOMIAL TERMS
end
for i=1:k
y_hat(i)=sum(T(i,:));
end
Sr=sum((Y-y_hat).^2);
Se=sqrt(Sr/(k-(n+1))); %STANDARD ERROR-Se

%COEFFICIENT OF DETERMINATION
Cd=(St-Sr)/St %COEFFICIENT OF DETERMINATION (r^2)
Cd = 1.0000
fprintf(‘For n=%d. Coefficient of Determination=%0.5fn’,n,Cd)

function [dydt] = Projectfunction(t,x,y)
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

x1 = x(1); %x1=Xf
x2 = x(2); %x2=Xr
x3 = x(3); %x3=Xb
x4 = x(4); %x4=theta

y1 = y(1); %y1=Xf’=dXf/dt
y2 = y(2); %y2=Xr’=dXr/dt
y3 = y(3); %y3=Xb’=dXb/dt
y4 = y(4); %y4=theta’=dtheta/dt

dydt1 = (Ksf*((x3-(L1*x4)))-x1)+Bsf*((y3-(L1*y4)-y1)-(Kf*x1))/Mf;
dydt2 = (Ksr*((x3+(L2*x4))-x2)+Bsr*((y3+(L2*y4))-y2)-(Kr*x2))/Mr;
dydt3 = (Ksf*((x3-(L1*x4))-x1)+Bsf*((y3-(L1*y4))-y1)+Ksr*((x3+(L2*x4))-x2)+Bsr*((y3+(L2*y4))-y2)-(10*exp(-(5*t))))/Mb;
dydt4 = ((-(Ksf*(x1-(L1*x4))*L1)-(Bsf*(y1-(L1*y4))*L1)+(Ksr*(x2+(L2*x4))*L2)+(Bsr*(y2+(L2*y4))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
[dydt] = [dydt1; dydt2; dydt3; dydt4];
end

% 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]};I get this error when trying to solve these 4 second order differential equations using a program that solves them using Heun’s method and then Gaussian elimination to get a polynomial equation the represents the system. Previsouly posted for help, but the soluitions have symbolics and I need a singular soluition that represents all 4 equations. Does any one have a suggestion on how I can slove this? Is using the "function" with all 4 equations the best way to do this, or is there a nicer method for doing this?
clear,close,clc