I Can't get this code to work, can anyone help

I'm working on a school project and we have to solve a quarter car physical model, and my code keeps giving me the same error for ode45

These are what I've created

%% P3
clc
for bs=[1000 1200 2000]
    odefun=@(t,xg)p3SPJ19(t,xg,bs);
    tspan=[0 10];
    xg=[0.01 0 0 0];
    [t,xg]=ode45(odefun,tspan,xg);
    hold on
    subplot(2,1,1)
    plot(t,xg(:,1));
    subplot(2,1,2)
    plot(t,xg(:,2));
end
hold off
legend('bs=1000')

function dxg = p3SPJ19 (t,xg,bs)
    kt=135000;
    ks=5700;
    bt=1400;
    bs=[1000 1200 2000];
    mc=466.5;
    mw=49.8;
    syms xw xw(2) xw(3) xc xc(2) xg xg(2)
    xg=((mw*xw3+bs*(xw2-xc2)+ks*(xw-xc)-bt*(xg2-xw2))/kt)+xw;
    dxg =zeroes(4,1);
    dxg(1)=xg(3);
    xg(3)=xc*(xw - (bs*(xc2 - xw2) - mw*xw3 + bt*(xg2 - xw2))/kt) - (ks*(xc - xw)^2)/(2*kt);
    dxg(2)=xg(4);
    xg(4)=- (ks*xc^3)/(6*kt) + (xw/2 - (bs*(xc(2) - xw(2)) - mw*xw(3) + bt*(xg(2) - xw(2)))/(2*kt) + (ks*xw)/(2*kt))*xc^2 - (ks*xw^2*xc)/(2*kt);
    dxg(3)=- (ks*xc^4)/(24*kt) + (xw/6 - (bs*(xc(2) - xw(2)) - mw*xw(3) + bt*(xg(2) - xw(2)))/(6*kt) + (ks*xw)/(6*kt))*xc^3 - (ks*xw^2*xc^2)/(4*kt);
    dxg(4)=- (ks*xc^5)/(120*kt) + (xw/24 - (bs*(xc(2) - xw(2)) - mw*xw(3) + bt*(xg(2) - xw(2)))/(24*kt) + (ks*xw)/(24*kt))*xc^4 - (ks*xw^2*xc^3)/(12*kt);
end

I think I may have messed up my function

I'm clueless on the Simulink portion as well

-----

I started over on the function, is this any better?

function xg = p3SPJ19 (t,xg,bs)
kt=135000;
  ks=5700;
  bt=1400;
  bs=[1000 1200 2000];
  mc=466.5;
  mw=49.8;
xc=-((mc*xw(2)+bs*(xw(1)-xc(1)))/(kt))+(xw);
xc(1)=diff(xc);
xg=zeros(2,1);
xg(0)=((mw*xw(3)+bs*(xw(2)-xc(2))+ks*(xw-xc)-bt*(xg(2)-xw(2)))/kt)+xw;
xg(1)=xc*(xw - (bs*(xc(2) - xw(2)) - mw*xw(3) + bt*(xg(2) - xw(2)))/kt) - (ks*(xc - xw)^2)/(2*kt);
xg(2)=- (ks*xc^3)/(6*kt) + (xw/2 - (bs*(xc(2) - xw(2)) - mw*xw(3) + bt*(xg(2) - xw(2)))/((2)*kt) + (ks*xw)/(2*kt))*xc^2 - (ks*xw^2*xc)/(2*kt);

[James Tursa]

ode45 is a numeric solver. Why are you using symbolic variables inside a numeric solver? You need to calculate dxg based on the inputs t, xg, and bs.

[Walter Roberson]

See the documentation for odeFunction to turn an ode into a numeric system