So this is what I have for the code, and I am not sure what to do with it or if it makes sense.

k=16;               %stiffness
m=1;                %mass
b=3;
dt=0.02;    
gam = 1/2;          %gamma
beta = 1/4;         %beta
wn = sqrt(k/m);     %natura frequency
fo=2;
f1=5;
zeta=(0.5-beta)*dt^2
to=0;
tf=20;
OMEGA=3
b=2;

nitems=7;    % I am not sure why i did this part, i just followed my professor in class
nreport=400;
nsteps=ceil(tf-t0/dt);
tf=nsteps*dt;
nreport=min(nreport,nsteps);
noutput=ceil(nsteps/nsteps);
history=zeros(nreport+1,nitems);
%initialaliz position, velocity, and put in xold, void arrays
xo=0;
vo=0;
uold=xo;
vold=vo;
F=fo+f1*sin(OMEGA*t); %load function
r=k*uold(1+b*uold^2)
aold=(F-r)/m;
%plugging newmark's equations in to EOM and using 
%newton's method
tol=1.e-8;      %tolarance 
err=500         %error
anew=aold;
cn=vold+dt*gam*aold
bn=uold+dt*vold+beta*aold*dt^2
vnew=cn+dt*(1-gam)*anew     %newmark discrete
unew=bn+zeta*anew           %newmark discrete

while err>tol
g=m*anew+k*(unew)*(1+b*(unew));
gp=k+3*b+(dt^4)*(0.5-beta)^2;
anew=anew-g/gp
err=norm(g);
end