Quick questions on mechanical vibrations

a)Use the same model as in first case but for c=0.2 and let the driving frequency be constant ωdr=4 with zero initial conditions. Plot the motion of the system and describe the motion.

b) Let the force be accelerating such that (F sin(5t^2 /1000)). Plot the response (for c=0.2) for each mass as function of angular frequency. Explain the result

I have already solved for part a, however for part b, it requires me to plot the response for each mass as a function of angular frequency, does any one have any idea on how do I go about it? I have already changed the force F to the required parameters.

Main code:

% Simulation of forced 2DOF system with ode45
clear all
clc
global A B F Z I wdr
c=0.2;
m=1;
k=100;
c=0.2;
F0=10;
t0=0;
t1=100;
M=[m 0 0 ;0 m 0 ; 0 0 m];
K=[2*k -k 0;-k 2*k -k;0 -k k];
C=[2*c -c 0;-c 2*c -c;0 -c c];
for t=t0:t1
    f=[0 0 F0*sin(5*t^2/1000)]';
end
% Force vector

wdr=4; % Driving frequency
A=M\K;
B=M\C;
F=M\f;
Z=zeros(3); % Zero matrix
I=eye(3); % Diagonal matrix
x0=[0 0 0 0 0 0]'; % Initial conditions
[t,x]=ode45('Fun1',[t0 t1],x0);
dx=[t,x]

plot(t,x(:,1),t,x(:,2));

Func 1 code:

function dx=Fun2(t,x)
global A B F Z I wdr
dx=zeros(6,1);

dx=[Z I ; -A -B] *x+[[0 0 0]' ; F]*sin(wdr*t);

end

Answer 1

[KSSV]

Check the below code......you should run a loop for each mass. It will give output as R. It has the required results of each mass.

function R =  myfunction()

global A B F Z I wdr
c=0.2;
mass=1:5;  % masses 
% loop for each mass 
R = cell(length(mass),1) ;  %Result 
for i = 1:length(mass)
    m = mass(i) ;
    k=100;
    c=0.2;
    F0=10;
    t0=0;
    t1=100;
    M=[m 0 0 ;0 m 0 ; 0 0 m];
    K=[2*k -k 0;-k 2*k -k;0 -k k];
    C=[2*c -c 0;-c 2*c -c;0 -c c];
    for t=t0:t1
        f=[0 0 F0*sin(5*t^2/1000)]';
    end
    % Force vector
    wdr=4; % Driving frequency
    A=M\K;
    B=M\C;
    F=M\f;
    Z=zeros(3); % Zero matrix
    I=eye(3); % Diagonal matrix
    x0=[0 0 0 0 0 0]'; % Initial conditions
    [t,x]=ode45(@Fun1,[t0 t1],x0);
    R{i}=[t,x] ;
    plot(t,x(:,1),t,x(:,2));
    title(sprintf('mass = %s',num2str(m)))
    drawnow
end
end
function dx=Fun1(t,x)
global A B F Z I wdr
dx=zeros(6,1);
dx=[Z I ; -A -B] *x+[[0 0 0]' ; F]*sin(wdr*t);
end

[yung zhen lek]

Thanks for the help, greatly appreciated =)