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

% 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]};

%-------------------------------------------------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.2f\t y=%0.3f\t z=%0.3f\n',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=A\B                           %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.5f\n',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]};

[Torsten]

Starting from the line

%%%LEAST SQUARE METHOD-FINDS POLYNOMIAL FOR GIVEN DATA SET%%%%%

I have no idea of what the code is supposed to do.

[John D'Errico]

I see MANY issues here. You have n data points. Actually n+1, since 0 is one of your points. Then you want to use least squares to solve for a polynomial. where the polynomial has n+1 coefficients. As such, the polynomial will be an interpolating polynomial. Least squares is not at all relevant, since an interpolating poynomial will pass exactly through the data points.

Then you go through a great deal of confusion to use both Gaussian elimination, and also backslash?? JUST USE POLYFIT FOR GOD SAKES!

You actually generate the normal equations for a square system. Why square the condition number of a probem that for a 4th degree polynomial will already be difficult?

Computing standard errors is a waste of time for an interpolating polynomial.

I could go on and on of course.

[John D'Errico]

To be honest, I'd want to use a very different method. That is, start by assuming the solution can be represented by a series as an approximation. Now use minimum energy methods to find the series solution to the differential equations. The result will be a polynomial, exactly as you want. Another approach would be to use the chebfun toolbox to solve the problem. Again, you will have a "function" as the solution. Either approach would be greatly preferable in my eyes.

[Jon]

For this project I am supposed to use 2 different numerical methods to find the polynomial. Although using "POLYFIT" would be easier that doen not fit what I am supposed to do. That is why I am tryong to find out how to solve the 4 dependent second-order differential equations. I thought uning "Function" but I get the errors above. It has also been suggested by a classmate to use "mymatlabFunction" but they didn't know how I could use it. That is what I am trying to figure out which would work best within the scope of the assignment.

[Torsten]

    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;
    

This is not Heun's method.

It cannot happen that xn(i+1) and yn(i+1) appear on both sides of the updates because Heun's method is explicit.