trying to find the polynomial equation representing the given four second-order differential equations

For the following problem in the first image i get the error pictured in the third image when trying to solve the 4 second-order differential equations using numerical methods. I am trying to use Heun's method for solving second-order differential equations and from here getting the x and y values corresponding to the function. Then taking these values and putting them into the least-square method to find the polynomial equation that best represents the system of equations. The third image is of the free-body diagram representing this system, then follows the code, if anyone can see where I need to adjust the code to run correctly this would be much appreciated.

MATLAB CODE:

% 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

%Initial Conditions x=(0)=0 dXf/dt(0)=y(0)=0  ;
%t from 0 to 1 h=dx=0.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

syms Xf(t) Xr(t) Xb(t) theta(t) fa(t)  Y

Xf_1=diff(Xf,t);
Xf_2=diff(Xf,t,2);

Xr_1=diff(Xr,t);
Xr_2=diff(Xr,t,2);

Xb_1=diff(Xb,t);
Xb_2=diff(Xb,t,2);

theta_1=diff(theta,t);
theta_2=diff(theta,t,2);

%%%% MfXf"=Ksf((Xb-(L1*theta))-Xf)+Bsf((Xb'-(L1*theta'))-Xf')-(Kf*Xf);

%anything in [] below indicates a deritive
% Mf*[(d^2Xf)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)])-[(dXf)/(dt)])-(Kf*Xf);

Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0;
f2_YF=isolate(Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0,Xf_2);

f2_Yf=subs(f2_YF,[Xf_2 Xf_1],[Xf_1 Y])

f2_Yf = 

f2_Yf=rhs(f2_Yf)

f2_Yf = 

%dXf/dt=y_f-->ASSUMED
%dXf/dt=f1_f(x,y,t)=y_f
%dYf/dt=f2_f(x,y,t)=f2_Yf

%%%% MrXr"=Ksr((Xb+(L2*theta))-Xr)+Bsr((Xb'+(L2*theta'))-Xr')-(Kr*Xr);

%anything in [] below indicates a deritive
% Mf*[(d^2Xr)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsr*(([(dXb)/(dt)]-(L2*[(dtheta)/(dt)])-[(dXr)/(dt)])-(Kf*Xf);
Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0;

f2_YR=isolate(Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0,Xr_2);

f2_Yr=subs(f2_YR,[Xr_2 Xr_1],[Xr_1 Y])

f2_Yr = 

f2_Yr=rhs(f2_Yr)

f2_Yr = 

%dXr/dt=y_r-->ASSUMED
%dXr/dt=f1_r(x,y,t)=y_r
%dYr/dt=f2_r(x,y,t)=f2_Yr

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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
% Mb*[(d^2Xb)/(dt^2)]=Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)]))-[(dXf)/(dt)])+Ksr*([Xb+(L2*theta)]-Xr)+Bsr( (([(dXb)/(dt)]+(L2*[(dtheta)/(dt)]))-[(dXr)/(dt)]))+fa(t);
Mb*Xb_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1))-Xf_1)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1))==(10*exp(-(5*t))));
f2_YB=Xb_2==(-Ksf*((Xb-(L1*theta))-Xf)-Bsf*((Xb_1-(L1*theta_1))-Xf_1)-Ksr*((Xb+(L2*theta))-Xr)-Bsr*((Xb_1+(L2*theta_1))-Xr_1)+(10*exp(-(5*t))))/Mb;

f2_Yb=subs(f2_YB,[Xb_2 Xb_1],[Xb_1 Y])

f2_Yb(t) = 

f2_Yb=rhs(f2_Yb)

f2_Yb(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_b(x,y,t)=y_b
%dYb/dt=f2_b(x,y,t)=f2_Yb



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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
%Ic*[(d^2theta)/(dt^2)]=(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*([dXf/dt]-(L1*[(dtheta)/(dt)]))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*([(dXr)/(dt)]+(L2*[(dtheta)/(dt)]))*L2)+(fa(t)*L3));

theta_2==((-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
Ic*theta_2-(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2))==(10*exp(-(5*t)))*L3;


f2_YTheta=theta_2==((-Ksf*(Xf-(Xb-(L1*theta)))*L1) - (Bsf*(Xf_1-(Xb_1-(L1*theta_1)))*L1) + (Ksr*(Xr-(Xb+(L2*theta)))*L2) + (Bsr*(Xr_1-(Xb_1+(L2*theta_1)))*L2) + ((10*exp(-(5*t)))*L3))/Ic;

f2_Ytheta=subs(f2_YTheta,[theta_2 theta_1],[theta_1 Y])

f2_Ytheta(t) = 

f2_Ytheta=rhs(f2_Ytheta)

f2_Ytheta(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_theta(x,y,t)=y_theta
%dYb/dt=f2_theta(x,y,t)=f2_Ytheta

syms x  y_f  y_r  y_b  y_theta  t

%==INPUT SECTION==%

fx_f=@(x,y,t)y_f;

fy_f=@(x,y,t)f2_Yf;

fx_r=@(x,y,t)y_r;
fy_r=@(x,y,t)f2_Yr;

fx_b=@(x,y,t)y_b;
fy_b=@(x,y,t)f2_Yb;

fx_theta=@(x,y,t)y_theta;
fy_theta=@(x,y,t)f2_Ytheta;


%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn_f(1)=x0;
yn_f=y0;
xn_r(1)=x0;
yn_r=y0;
xn_b(1)=x0;
yn_b=y0;
xn_theta(1)=x0;
yn_theta=y0;


for i=1:length(tn)
     %==EULER'S METHOD
    xn_f(i+1)=xn_f(i)+fx_f(xn_f(i),yn_f(i),tn(i))*h;
    yn_f(i+1)=yn_f(i)+fy_f(xn_f(i),yn_f(i),tn(i))*h;

    xn_r(i+1)=xn_r(i)+fx_r(xn_r(i),yn_r(i),tn(i))*h;
    yn_r(i+1)=yn_r(i)+fy_r(xn_r(i),yn_r(i),tn(i))*h;

    xn_b(i+1)=xn_b(i)+fx_b(xn_b(i),yn_b(i),tn(i))*h;
    yn_b(i+1)=yn_b(i)+fy_b(xn_b(i),yn_b(i),tn(i))*h;

    xn_theta(i+1)=xn_theta(i)+fx_theta(xn_theta(i),yn_theta(i),tn(i))*h;
    yn_theta(i+1)=yn_theta(i)+fy_theta(xn_theta(i),yn_theta(i),tn(i))*h;

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

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

    xn_f(i+1)=xn_f(i)+0.5*(fx_f(xn_f(i),yn_f(i),tn(i))+fx_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;
    yn_f(i+1)=yn_f(i)+0.5*(fy_f(xn_f(i),yn_f(i),tn(i))+fy_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;

    xn_r(i+1)=xn_r(i)+0.5*(fx_r(xn_r(i),yn_r(i),tn(i))+fx_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;
    yn_r(i+1)=yn_r(i)+0.5*(fy_r(xn_r(i),yn_r(i),tn(i))+fy_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;

    xn_b(i+1)=xn_b(i)+0.5*(fx_b(xn_b(i),yn_b(i),tn(i))+fx_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;
    yn_b(i+1)=yn_b(i)+0.5*(fy_b(xn_b(i),yn_b(i),tn(i))+fy_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;

    xn_theta(i+1)=xn_theta(i)+0.5*(fx_theta(xn_theta(i),yn_theta(i),tn(i))+fx_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;
    yn_theta(i+1)=yn_theta(i)+0.5*(fy_theta(xn_theta(i),yn_theta(i),tn(i))+fy_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;


    fprintf('t=%0.2f\t x_f=%0.3f\t y_f=%0.3f\t x_r=%0.3f\t y_r=%0.3f\t x_b=%0.3f\t y_b=%0.3f\t x_theta=%0.3f\t  y_theta=%0.3f\n',tn(i),xn_f(i),yn_f(i),xn_r(i),yn_r(i),xn_b(i),yn_b(i),xn_theta(i),yn_theta(i))
end

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using mupadengine/feval2char
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

%input
x=['x values from Heuns method'];
y=['y values from Heuns method'];
n=4;                   %order of the polynomial you want to find 

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

%INPUT SECTION
x=x;                   %x-cordinates from data-input; independent vairiables   
y=y;                   %y-cordinates from data-output; dependent vairiables



%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)
%========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)
fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)
%EQUATION FOR THE POLYNOMIAL
syms x y

for i=0:n
    for n=1 
        a0=a(:,1);a1=a(:,2);a2=a(:,3);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2
    end
    
    for n=2 
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3
    end

    for n=3 
        a0=a(:,1);a1=a(:,2);a2=a(:,2);a3=a(:,4);a4=a(:,5);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4
    end

    for n=4 
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5
    end

    for n=5 
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7)
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6
    end

    for n=6 
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7
    end

    for n=7
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8
    end

    for n=8
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9
    end

    for n=9
        a0=a(:,1);a1=a(:,2);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);a10=a(:,11);
        sympref("FloatingPointOutput",true);
        y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9+a10*x.^10
    end   
end

[Walter Roberson]

fx_f=@(x,y,t)y_f;

That code declares that fx_t is to be the handle to an anonymous function. The anonymous function accepts up to three inputs. The three inputs are ignored, and the constant output y_f is to be returned.

Except... you have not defined any y_f

We see a hint in

%dXf/dt=f1_f(x,y,t)=y_f

but there is no f1_f defined either.

[Sam Chak]

Since 'Xf(t)' is declared as a symbolic variable and the hint tells us that it is a rate of change of 'Xf', perhaps @Jon intends to write the code as follows:

y_f = diff(Xf);

[Jon]

For the method that I am supposed to use he has us changing the first order into a symbolic y. I did forget to add the symbolic y_f y_r y_b and y_theta into the syms command line. the code after finding f2_Ytheta should appear as follows:

% 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

%Initial Conditions x=(0)=0 dXf/dt(0)=y(0)=0  ;
%t from 0 to 1 h=dx=0.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

syms Xf(t) Xr(t) Xb(t) theta(t) fa(t)  Y

Xf_1=diff(Xf,t);
Xf_2=diff(Xf,t,2);

Xr_1=diff(Xr,t);
Xr_2=diff(Xr,t,2);

Xb_1=diff(Xb,t);
Xb_2=diff(Xb,t,2);

theta_1=diff(theta,t);
theta_2=diff(theta,t,2);

%%%% MfXf"=Ksf((Xb-(L1*theta))-Xf)+Bsf((Xb'-(L1*theta'))-Xf')-(Kf*Xf);

%anything in [] below indicates a deritive
% Mf*[(d^2Xf)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)])-[(dXf)/(dt)])-(Kf*Xf);

Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0;
f2_YF=isolate(Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0,Xf_2);

f2_Yf=subs(f2_YF,[Xf_2 Xf_1],[Xf_1 Y])

f2_Yf = 

f2_Yf=rhs(f2_Yf)

f2_Yf = 

%dXf/dt=y_f-->ASSUMED
%dXf/dt=f1_f(x,y,t)=y_f
%dYf/dt=f2_f(x,y,t)=f2_Yf

%%%% MrXr"=Ksr((Xb+(L2*theta))-Xr)+Bsr((Xb'+(L2*theta'))-Xr')-(Kr*Xr);

%anything in [] below indicates a deritive
% Mf*[(d^2Xr)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsr*(([(dXb)/(dt)]-(L2*[(dtheta)/(dt)])-[(dXr)/(dt)])-(Kf*Xf);
Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0;

f2_YR=isolate(Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0,Xr_2);

f2_Yr=subs(f2_YR,[Xr_2 Xr_1],[Xr_1 Y])

f2_Yr = 

f2_Yr=rhs(f2_Yr)

f2_Yr = 

%dXr/dt=y_r-->ASSUMED
%dXr/dt=f1_r(x,y,t)=y_r
%dYr/dt=f2_r(x,y,t)=f2_Yr

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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
% Mb*[(d^2Xb)/(dt^2)]=Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)]))-[(dXf)/(dt)])+Ksr*([Xb+(L2*theta)]-Xr)+Bsr( (([(dXb)/(dt)]+(L2*[(dtheta)/(dt)]))-[(dXr)/(dt)]))+fa(t);
Mb*Xb_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1))-Xf_1)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1))==(10*exp(-(5*t))));
f2_YB=Xb_2==(-Ksf*((Xb-(L1*theta))-Xf)-Bsf*((Xb_1-(L1*theta_1))-Xf_1)-Ksr*((Xb+(L2*theta))-Xr)-Bsr*((Xb_1+(L2*theta_1))-Xr_1)+(10*exp(-(5*t))))/Mb;

f2_Yb=subs(f2_YB,[Xb_2 Xb_1],[Xb_1 Y])

f2_Yb(t) = 

f2_Yb=rhs(f2_Yb)

f2_Yb(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_b(x,y,t)=y_b
%dYb/dt=f2_b(x,y,t)=f2_Yb



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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
%Ic*[(d^2theta)/(dt^2)]=(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*([dXf/dt]-(L1*[(dtheta)/(dt)]))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*([(dXr)/(dt)]+(L2*[(dtheta)/(dt)]))*L2)+(fa(t)*L3));

theta_2==((-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
Ic*theta_2-(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2))==(10*exp(-(5*t)))*L3;


f2_YTheta=theta_2==((-Ksf*(Xf-(Xb-(L1*theta)))*L1) - (Bsf*(Xf_1-(Xb_1-(L1*theta_1)))*L1) + (Ksr*(Xr-(Xb+(L2*theta)))*L2) + (Bsr*(Xr_1-(Xb_1+(L2*theta_1)))*L2) + ((10*exp(-(5*t)))*L3))/Ic;

f2_Ytheta=subs(f2_YTheta,[theta_2 theta_1],[theta_1 Y])

f2_Ytheta(t) = 

f2_Ytheta=rhs(f2_Ytheta)

f2_Ytheta(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_theta(x,y,t)=y_theta
%dYb/dt=f2_theta(x,y,t)=f2_Ytheta

syms x  y_f  y_r  y_b  y_theta  t

%==INPUT SECTION==%

fx_f=@(x,y,t)y_f;

fy_f=@(x,y,t)f2_Yf;

fx_r=@(x,y,t)y_r;
fy_r=@(x,y,t)f2_Yr;

fx_b=@(x,y,t)y_b;
fy_b=@(x,y,t)f2_Yb;

fx_theta=@(x,y,t)y_theta;
fy_theta=@(x,y,t)f2_Ytheta;

%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn_f(1)=x0;
yn_f=y0;
xn_r(1)=x0;
yn_r=y0;
xn_b(1)=x0;
yn_b=y0;
xn_theta(1)=x0;
yn_theta=y0;


for i=1:length(tn)
    %==EULER'S METHOD
    xn_f(i+1)=xn_f(i)+fx_f(xn_f(i),yn_f(i),tn(i))*h;
    yn_f(i+1)=yn_f(i)+fy_f(xn_f(i),yn_f(i),tn(i))*h;
    
    xn_r(i+1)=xn_r(i)+fx_r(xn_r(i),yn_r(i),tn(i))*h;
    yn_r(i+1)=yn_r(i)+fy_r(xn_r(i),yn_r(i),tn(i))*h;
    
    xn_b(i+1)=xn_b(i)+fx_b(xn_b(i),yn_b(i),tn(i))*h;
    yn_b(i+1)=yn_b(i)+fy_b(xn_b(i),yn_b(i),tn(i))*h;
    
    xn_theta(i+1)=xn_theta(i)+fx_theta(xn_theta(i),yn_theta(i),tn(i))*h;
    yn_theta(i+1)=yn_theta(i)+fy_theta(xn_theta(i),yn_theta(i),tn(i))*h;
    
    %==NEXT 3 LINES ARE FOR HEUN'S METHOD
    
    tn(i+1)=tn(i)+h;
    
    xn_f(i+1)=xn_f(i)+0.5*(fx_f(xn_f(i),yn_f(i),tn(i))+fx_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;
    yn_f(i+1)=yn_f(i)+0.5*(fy_f(xn_f(i),yn_f(i),tn(i))+fy_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;
    
    xn_r(i+1)=xn_r(i)+0.5*(fx_r(xn_r(i),yn_r(i),tn(i))+fx_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;
    yn_r(i+1)=yn_r(i)+0.5*(fy_r(xn_r(i),yn_r(i),tn(i))+fy_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;
    
    xn_b(i+1)=xn_b(i)+0.5*(fx_b(xn_b(i),yn_b(i),tn(i))+fx_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;
    yn_b(i+1)=yn_b(i)+0.5*(fy_b(xn_b(i),yn_b(i),tn(i))+fy_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;
    
    xn_theta(i+1)=xn_theta(i)+0.5*(fx_theta(xn_theta(i),yn_theta(i),tn(i))+fx_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;
    yn_theta(i+1)=yn_theta(i)+0.5*(fy_theta(xn_theta(i),yn_theta(i),tn(i))+fy_theta(xn-theta(i+1),yn_theta(i+1),tn(i+1)))*h;
    
    
    fprintf('t=%0.2f\t x_f=%0.3f\t y_f=%0.3f\t x_r=%0.3f\t y_r=%0.3f\t x_b=%0.3f\t y_b=%0.3f\t x_theta=%0.3f\t  y_theta=%0.3f\n',tn(i),xn_f(i),yn_f(i),xn_r(i),yn_r(i),xn_b(i),yn_b(i),xn_theta(i),yn_theta(i))
end

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using mupadengine/feval2char
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

%input
x=['x values from Heuns method'];
y=['y values from Heuns method'];
n=4;                   %order of the polynomial you want to find 

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

%INPUT SECTION
x=x;                   %x-cordinates from data-input; independent vairiables   
y=y;                   %y-cordinates from data-output; dependent vairiables



%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)
%========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)
fprintf('For n=%d. Coefficient of Determination=%0.5f\n',n,Cd)
%EQUATION FOR THE POLYNOMIAL
syms x y
%for n=2 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2
%for n=3 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3
%for n=4 activate lines below
a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4
%for n=5 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5
%for n=6 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7)
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6
%for n=7 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7
%for n=8 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8
%for n=9 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9

but doing that I get the following error. Any suggestions

[Sam Chak]

Hi @Jon, Would you edit the post and click the Indentation icon to paste the MATLAB code in proper format?

Clicking the icon to toggle the "grey field" that interprets MATLAB code when the Play (Run) icon is clicked.

%% Paste your code in this grey field

[Torsten]

This part of your code will not work:

syms x  y_f  y_r  y_b  y_theta  t
%==INPUT SECTION==%
fx_f=@(x,y,t)y_f;
fy_f=@(x,y,t)f2_Yf;
fx_r=@(x,y,t)y_r;
fy_r=@(x,y,t)f2_Yr;
fx_b=@(x,y,t)y_b;
fy_b=@(x,y,t)f2_Yb;
fx_theta=@(x,y,t)y_theta;
fy_theta=@(x,y,t)f2_Ytheta;

Take a look here about the relevant commands to convert symbolic differential equations to a form that is suited for a numerical solver:

https://blogs.mathworks.com/loren/2013/06/10/from-symbolic-differential-equations-to-their-numeric-solution/

[Sam Chak]

Hi @Jon,

The image of the Mass-Spring-Damper Lever diagram is quite blurry. If possible, please upload a high-quality image.

Since you already have the system of differential equations, why not straightaway code them in the Linear State-Space model?

where is a state matrix and is a input matrix.

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

[Jon]

I did that for another class, for this class we are supposed to use numerical methods to come up with the polynomial that represents this model.

[Walter Roberson]

% 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

%Initial Conditions x=(0)=0 dXf/dt(0)=y(0)=0  ;
%t from 0 to 1 h=dx=0.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

syms Xf(t) Xr(t) Xb(t) theta(t) fa(t)  Y

Xf_1=diff(Xf,t);
Xf_2=diff(Xf,t,2);

Xr_1=diff(Xr,t);
Xr_2=diff(Xr,t,2);

Xb_1=diff(Xb,t);
Xb_2=diff(Xb,t,2);

theta_1=diff(theta,t);
theta_2=diff(theta,t,2);

%%%% MfXf"=Ksf((Xb-(L1*theta))-Xf)+Bsf((Xb'-(L1*theta'))-Xf')-(Kf*Xf);

%anything in [] below indicates a deritive
% Mf*[(d^2Xf)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)])-[(dXf)/(dt)])-(Kf*Xf);

Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0;
f2_YF=isolate(Mf*Xf_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1)-Xf_1)-(Kf*Xf))))==0,Xf_2);

f2_Yf=subs(f2_YF,[Xf_2 Xf_1],[Xf_1 Y])

f2_Yf = 

f2_Yf=rhs(f2_Yf)

f2_Yf = 

%dXf/dt=y_f-->ASSUMED
%dXf/dt=f1_f(x,y,t)=y_f
%dYf/dt=f2_f(x,y,t)=f2_Yf

%%%% MrXr"=Ksr((Xb+(L2*theta))-Xr)+Bsr((Xb'+(L2*theta'))-Xr')-(Kr*Xr);

%anything in [] below indicates a deritive
% Mf*[(d^2Xr)/(dt^2)] = Ksf*((Xb-(L1*theta))-Xf)+Bsr*(([(dXb)/(dt)]-(L2*[(dtheta)/(dt)])-[(dXr)/(dt)])-(Kf*Xf);
Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0;

f2_YR=isolate(Mr*Xr_2-((Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1)-(Kr*Xr)))==0,Xr_2);

f2_Yr=subs(f2_YR,[Xr_2 Xr_1],[Xr_1 Y])

f2_Yr = 

f2_Yr=rhs(f2_Yr)

f2_Yr = 

%dXr/dt=y_r-->ASSUMED
%dXr/dt=f1_r(x,y,t)=y_r
%dYr/dt=f2_r(x,y,t)=f2_Yr

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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
% Mb*[(d^2Xb)/(dt^2)]=Ksf*((Xb-(L1*theta))-Xf)+Bsf*(([(dXb)/(dt)]-(L1*[(dtheta)/(dt)]))-[(dXf)/(dt)])+Ksr*([Xb+(L2*theta)]-Xr)+Bsr( (([(dXb)/(dt)]+(L2*[(dtheta)/(dt)]))-[(dXr)/(dt)]))+fa(t);
Mb*Xb_2-((Ksf*((Xb-(L1*theta))-Xf)+Bsf*((Xb_1-(L1*theta_1))-Xf_1)+Ksr*((Xb+(L2*theta))-Xr)+Bsr*((Xb_1+(L2*theta_1))-Xr_1))==(10*exp(-(5*t))));
f2_YB=Xb_2==(-Ksf*((Xb-(L1*theta))-Xf)-Bsf*((Xb_1-(L1*theta_1))-Xf_1)-Ksr*((Xb+(L2*theta))-Xr)-Bsr*((Xb_1+(L2*theta_1))-Xr_1)+(10*exp(-(5*t))))/Mb;

f2_Yb=subs(f2_YB,[Xb_2 Xb_1],[Xb_1 Y])

f2_Yb(t) = 

f2_Yb=rhs(f2_Yb)

f2_Yb(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_b(x,y,t)=y_b
%dYb/dt=f2_b(x,y,t)=f2_Yb



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

%anything in [] below indicates a deritive  fa(t) given as 10e^(-(5t))
%Ic*[(d^2theta)/(dt^2)]=(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*([dXf/dt]-(L1*[(dtheta)/(dt)]))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*([(dXr)/(dt)]+(L2*[(dtheta)/(dt)]))*L2)+(fa(t)*L3));

theta_2==((-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2)+((10*exp(-(5*t)))*L3)))/Ic;
Ic*theta_2-(-(Ksf*(Xf-(L1*theta))*L1)-(Bsf*(Xf_1-(L1*theta_1))*L1)+(Ksr*(Xr+(L2*theta))*L2)+(Bsr*(Xr_1+(L2*theta_1))*L2))==(10*exp(-(5*t)))*L3;


f2_YTheta=theta_2==((-Ksf*(Xf-(Xb-(L1*theta)))*L1) - (Bsf*(Xf_1-(Xb_1-(L1*theta_1)))*L1) + (Ksr*(Xr-(Xb+(L2*theta)))*L2) + (Bsr*(Xr_1-(Xb_1+(L2*theta_1)))*L2) + ((10*exp(-(5*t)))*L3))/Ic;

f2_Ytheta=subs(f2_YTheta,[theta_2 theta_1],[theta_1 Y])

f2_Ytheta(t) = 

f2_Ytheta=rhs(f2_Ytheta)

f2_Ytheta(t) = 

%dXb/dt=y-->ASSUMED
%dXb/dt=f1_theta(x,y,t)=y_theta
%dYb/dt=f2_theta(x,y,t)=f2_Ytheta

syms x  y_f  y_r  y_b  y_theta  t

%==INPUT SECTION==%

fx_f=@(x,y,t)y_f;

fy_f=@(x,y,t)f2_Yf;

fx_r=@(x,y,t)y_r;
fy_r=@(x,y,t)f2_Yr;

fx_b=@(x,y,t)y_b;
fy_b=@(x,y,t)f2_Yb;

fx_theta=@(x,y,t)y_theta;
fy_theta=@(x,y,t)f2_Ytheta;

%==CALCULATIONS SECTION==%
tn=t0:h:tm;
xn_f(1) = sym(x0);
yn_f = sym(y0);
xn_r(1) = sym(x0);
yn_r = sym(y0);
xn_b(1) = sym(x0);
yn_b = sym(y0);
xn_theta(1) = sym(x0);
yn_theta = sym(y0);


for i=1:length(tn)
    %==EULER'S METHOD
    xn_f(i+1)=xn_f(i)+fx_f(xn_f(i),yn_f(i),tn(i))*h;
    yn_f(i+1)=yn_f(i)+fy_f(xn_f(i),yn_f(i),tn(i))*h;
    
    xn_r(i+1)=xn_r(i)+fx_r(xn_r(i),yn_r(i),tn(i))*h;
    yn_r(i+1)=yn_r(i)+fy_r(xn_r(i),yn_r(i),tn(i))*h;
    
    xn_b(i+1)=xn_b(i)+fx_b(xn_b(i),yn_b(i),tn(i))*h;
    yn_b(i+1)=yn_b(i)+fy_b(xn_b(i),yn_b(i),tn(i))*h;
    
    xn_theta(i+1)=xn_theta(i)+fx_theta(xn_theta(i),yn_theta(i),tn(i))*h;
    yn_theta(i+1)=yn_theta(i)+fy_theta(xn_theta(i),yn_theta(i),tn(i))*h;
    
    %==NEXT 3 LINES ARE FOR HEUN'S METHOD
    
    tn(i+1)=tn(i)+h;
    
    xn_f(i+1)=xn_f(i)+0.5*(fx_f(xn_f(i),yn_f(i),tn(i))+fx_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;
    yn_f(i+1)=yn_f(i)+0.5*(fy_f(xn_f(i),yn_f(i),tn(i))+fy_f(xn_f(i+1),yn_f(i+1),tn(i+1)))*h;
    
    xn_r(i+1)=xn_r(i)+0.5*(fx_r(xn_r(i),yn_r(i),tn(i))+fx_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;
    yn_r(i+1)=yn_r(i)+0.5*(fy_r(xn_r(i),yn_r(i),tn(i))+fy_r(xn_r(i+1),yn_r(i+1),tn(i+1)))*h;
    
    xn_b(i+1)=xn_b(i)+0.5*(fx_b(xn_b(i),yn_b(i),tn(i))+fx_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;
    yn_b(i+1)=yn_b(i)+0.5*(fy_b(xn_b(i),yn_b(i),tn(i))+fy_b(xn_b(i+1),yn_b(i+1),tn(i+1)))*h;
    
    xn_theta(i+1)=xn_theta(i)+0.5*(fx_theta(xn_theta(i),yn_theta(i),tn(i))+fx_theta(xn_theta(i+1),yn_theta(i+1),tn(i+1)))*h;
    yn_theta(i+1)=yn_theta(i)+0.5*(fy_theta(xn_theta(i),yn_theta(i),tn(i))+fy_theta(xn_theta(i+1),yn_theta(i+1),tn(i+1)))*h;
    
    
    fprintf('t=%s\t x_f=%s\t y_f=%s\t x_r=%s\t y_r=%s\t x_b=%s\t y_b=%s\t x_theta=%s\t  y_theta=%s\n',string(tn(i)),string(xn_f(i)),string(yn_f(i)),string(xn_r(i)),string(yn_r(i)),string(xn_b(i)),string(yn_b(i)),string(xn_theta(i)),string(yn_theta(i)))
end

t=0	 x_f=0	 y_f=0	 x_r=0	 y_r=0	 x_b=0	 y_b=0	 x_theta=0	  y_theta=0
t=5	 x_f=5*y_f	 y_f=(93750*Xb(t))/59 - (500*Y)/59 - (11593750*Xf(t))/59 - (271875*theta(t))/118 + (500*diff(Xb(t), t))/59 - (725*diff(theta(t), t))/59	 x_r=5*y_r	 y_r=(12574*Xb(t))/9 - (100*Y)/9 - (28756*Xr(t))/9 + (873893*theta(t))/450 + (100*diff(Xb(t), t))/9 + (139*diff(theta(t), t))/9	 x_b=5*y_b	 y_b=(5*exp(-5*t))/73 - (100*Y)/73 - (15662*Xb(t))/73 + (9375*Xf(t))/73 + (6287*Xr(t))/73 + (242741*theta(t))/3650 + (50*diff(Xf(t), t))/73 + (50*diff(Xr(t), t))/73 + (3*diff(theta(t), t))/73	 x_theta=5*y_theta	  y_theta=(67*exp(-5*t))/2712 - (20173*Y)/13560 + (242741*Xb(t))/6780 - (90625*Xf(t))/904 + (873893*Xr(t))/13560 - (159290251*theta(t))/678000 + (5*diff(Xb(t), t))/226 - (725*diff(Xf(t), t))/1356 + (695*diff(Xr(t), t))/1356
t=10	 x_f=10*y_f	 y_f=(187500*Xb(t))/59 - (1000*Y)/59 - (23187500*Xf(t))/59 - (271875*theta(t))/59 + (1000*diff(Xb(t), t))/59 - (1450*diff(theta(t), t))/59	 x_r=10*y_r	 y_r=(25148*Xb(t))/9 - (200*Y)/9 - (57512*Xr(t))/9 + (873893*theta(t))/225 + (200*diff(Xb(t), t))/9 + (278*diff(theta(t), t))/9	 x_b=10*y_b	 y_b=(10*exp(-5*t))/73 - (200*Y)/73 - (31324*Xb(t))/73 + (18750*Xf(t))/73 + (12574*Xr(t))/73 + (242741*theta(t))/1825 + (100*diff(Xf(t), t))/73 + (100*diff(Xr(t), t))/73 + (6*diff(theta(t), t))/73	 x_theta=10*y_theta	  y_theta=(67*exp(-5*t))/1356 - (20173*Y)/6780 + (242741*Xb(t))/3390 - (90625*Xf(t))/452 + (873893*Xr(t))/6780 - (159290251*theta(t))/339000 + (5*diff(Xb(t), t))/113 - (725*diff(Xf(t), t))/678 + (695*diff(Xr(t), t))/678
t=15	 x_f=15*y_f	 y_f=(281250*Xb(t))/59 - (1500*Y)/59 - (34781250*Xf(t))/59 - (815625*theta(t))/118 + (1500*diff(Xb(t), t))/59 - (2175*diff(theta(t), t))/59	 x_r=15*y_r	 y_r=(12574*Xb(t))/3 - (100*Y)/3 - (28756*Xr(t))/3 + (873893*theta(t))/150 + (100*diff(Xb(t), t))/3 + (139*diff(theta(t), t))/3	 x_b=15*y_b	 y_b=(15*exp(-5*t))/73 - (300*Y)/73 - (46986*Xb(t))/73 + (28125*Xf(t))/73 + (18861*Xr(t))/73 + (728223*theta(t))/3650 + (150*diff(Xf(t), t))/73 + (150*diff(Xr(t), t))/73 + (9*diff(theta(t), t))/73	 x_theta=15*y_theta	  y_theta=(67*exp(-5*t))/904 - (20173*Y)/4520 + (242741*Xb(t))/2260 - (271875*Xf(t))/904 + (873893*Xr(t))/4520 - (159290251*theta(t))/226000 + (15*diff(Xb(t), t))/226 - (725*diff(Xf(t), t))/452 + (695*diff(Xr(t), t))/452
t=20	 x_f=20*y_f	 y_f=(375000*Xb(t))/59 - (2000*Y)/59 - (46375000*Xf(t))/59 - (543750*theta(t))/59 + (2000*diff(Xb(t), t))/59 - (2900*diff(theta(t), t))/59	 x_r=20*y_r	 y_r=(50296*Xb(t))/9 - (400*Y)/9 - (115024*Xr(t))/9 + (1747786*theta(t))/225 + (400*diff(Xb(t), t))/9 + (556*diff(theta(t), t))/9	 x_b=20*y_b	 y_b=(20*exp(-5*t))/73 - (400*Y)/73 - (62648*Xb(t))/73 + (37500*Xf(t))/73 + (25148*Xr(t))/73 + (485482*theta(t))/1825 + (200*diff(Xf(t), t))/73 + (200*diff(Xr(t), t))/73 + (12*diff(theta(t), t))/73	 x_theta=20*y_theta	  y_theta=(67*exp(-5*t))/678 - (20173*Y)/3390 + (242741*Xb(t))/1695 - (90625*Xf(t))/226 + (873893*Xr(t))/3390 - (159290251*theta(t))/169500 + (10*diff(Xb(t), t))/113 - (725*diff(Xf(t), t))/339 + (695*diff(Xr(t), t))/339

%input
x=['x values from Heuns method'];
y=['y values from Heuns method'];
n=4;                   %order of the polynomial you want to find 

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

%INPUT SECTION
x=x;                   %x-cordinates from data-input; independent vairiables   
y=y;                   %y-cordinates from data-output; dependent vairiables



%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

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 5.426140e-22.

a1 = 5x1

27.3117 -0.7398 0.0371 -0.0003 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%%%%%=========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)

   27.3117   -0.7398    0.0371   -0.0003    0.0000

%========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)

For n=4. Coefficient of Determination=0.99998

%EQUATION FOR THE POLYNOMIAL
syms x y
%for n=2 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2
%for n=3 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3
%for n=4 activate lines below
a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);
sympref("FloatingPointOutput",true);
y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4

y = 

%for n=5 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5
%for n=6 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7)
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6
%for n=7 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7
%for n=8 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8
%for n=9 activate lines below
%a0=a(:,1);a1=a(:,1);a2=a(:,3);a3=a(:,4);a4=a(:,5);a5=a(:,6);a6=a(:,7);a7=a(:,8);a8=a(:,9);a9=a(:,10);
%sympref("FloatingPointOutput",true);
%y=a0+a1*x+a2*x.^2+a3*x.^3+a4*x.^4+a5*x.^5+a6*x.^6+a7*x.^7+a8*x.^8+a9*x.^9

[Jon]

thank you Walter, that was what I was trying to acheive.