%3b)
m=1000;  %mass must be >0
b=2*1000; %damper
k=40*1000; %spring

t0=0; % initial conditions
x0 = 0;
v0 = 0;

wn=sqrt(k/m);
zeta=[0.05,0.1,0.25,0.5,0.95];
U=[0:0.01:10];
opts=odeset('MaxStep',0.1);
maxamp=zeros(1,length(U));
maxforce=zeros(1,length(U));

for jj=1:length(zeta)

for ii=[1:size(U,2)]

A=[0,-1;wn^2,2*zeta(jj)*wn];
tspan = [t0 30]; %span in x variable
u0=[x0; v0];
[tarray,uarray] = ode45(@(t,u) -A*u+[0 ;(wn^2*r(U(ii)*t))], tspan, u0(:),opts); %ode45 ODE solver



maxamp(ii)=max(abs(uarray(:,1)));


% taking the variables


force=zeros(1,length(tarray));
for j=1:length(tarray)
  force(j)=-(k*(uarray(j,1)-r(U(ii)*tarray(j)))+(b*uarray(j,2)));


end
maxforce(ii)=max(abs(uarray(:,1)));


end
figure(jj+9)
plot(U,maxamp,'LineWidth',1.2,'DisplayName',"Resonse")
hold on
title("Max amplitude graph 3b)")
xlabel('U')
ylabel('max amplitude')
grid on

figure(jj+14)
plot(U,maxforce,'LineWidth',1.2,'DisplayName',"Resonse")
hold on
title("Max force graph 3b)")
xlabel('U')
ylabel('max force')
grid on

end
function rx=r(x)
Zr=0.1;
Lr=0.4;
xr=1;
if abs(x-xr)<Lr/2
    rx=-Zr;
else
    rx=0;
end
end