采用中心差分法计算其位移响应，并作出前10s内每个质点的时间位移曲线。m1=m2=1，m3=2，c1=c2=c3=1，k1=1,k2=5,k3=8,f1=0,f2=sin(w2*t),f3=0,w1=w2=w3=1,x1=x1'=x2=x2'=x3=x3'=0,matlab代码

根据题目中给出的物理模型，可以写出如下的运动方程：

m1 * x1'' + c1 * x1' + k1 * x1 - k2 * (x2 - x1) - f2 = 0
m2 * x2'' + c2 * x2' + k2 * (x2 - x1) - k3 * (x3 - x2) + f2 = 0
m3 * x3'' + c3 * x3' + k3 * (x3 - x2) = 0

其中，m1、m2、m3分别是三个质点的质量，c1、c2、c3分别是三个质点的阻尼系数，k1、k2、k3分别是三个质点之间的弹簧系数，f2是第二个质点外力的大小，w2是第二个质点外力的角频率，x1、x1'、x2、x2'、x3、x3'分别是三个质点的位移和速度。

为了计算三个质点的位移响应，需要将上述方程转化为一组一阶微分方程组。令

y1 = x1'
y2 = x2'
y3 = x3'

则有

x1' = y1
x2' = y2
x3' = y3
y1' = (-c1*y1 - k1*x1 + k2*(x2 - x1) + f2)/m1
y2' = (-c2*y2 - k2*(x2 - x1) + k3*(x3 - x2) - f2)/m2
y3' = (-c3*y3 - k3*(x3 - x2))/m3

下面是对应的MATLAB代码：

m1 = 1;
m2 = 1;
m3 = 2;
c1 = 1;
c2 = 1;
c3 = 1;
k1 = 1;
k2 = 5;
k3 = 8;
f1 = 0;
f2 = @(t) sin(t);
f3 = 0;
w1 = 1;
w2 = 1;
w3 = 1;
x10 = 0;
x20 = 0;
x30 = 0;
x10p = 0;
x20p = 0;
x30p = 0;

t = 0:0.01:10;
N = length(t);
x1 = zeros(1, N);
x2 = zeros(1, N);
x3 = zeros(1, N);

for i = 2:N
    y = [x10p, x20p, x30p, (-c1*x10p - k1*x10 + k2*(x20 - x10) + f2(t(i)))/m1, ...
         (-c2*x20p - k2*(x20 - x10) + k3*(x30 - x20) - f2(t(i)))/m2, ...
         (-c3*x30p - k3*(x30 - x20))/m3];
    k1 = dt * f(x(i-1), y);
    k2 = dt * f(x(i-1)+0.5*k1, y);
    k3 = dt * f(x(i-1)+0.5*k2, y);
    k4 = dt * f(x(i-1)+k3, y);
    x(i,:) = x(i-1,:) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    x10 = x(i,1);
    x20 = x(i,2);
    x30 = x(i,3);
    x10p = x(i,4);
    x20p = x(i,5);
    x30p = x(i,6);
end

figure;
plot(t, x(:,1), t, x(:,2), t, x(:,3));
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Displacement vs Time');
legend('x1', 'x2', 'x3');

figure;
plot(t, x(:,1), t, x(:,2), t, x(:,3), t, f2(t));
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Displacement vs Time');
legend('x1', 'x2', 'x3', 'f2');

其中，f是一个函数，用于计算微分方程组右侧的向量。