E = 2; rho = 1.0;
L = 10; a = 4; b = 0.1;
area = @(x) a - 2*b*x;

if area(L) < 0
    disp('A(L) < 0');
    exit();
end

n = 4; h = L/n;
nodal = h*[0:n];

for i=1:n
    Vs(i) = integral(area, nodal(i), nodal(i+1));
end
V = zeros(n+1,n+1);
for i=1:n
    V(i:i+1,i:i+1) = V(i:i+1,i:i+1) + Vs(i)*[1,-1;-1,1];
end
K = E/h/h*V;
K

% 解析解との比較
Kp = E/h*[
 a-b*h,    -a+b*h,         0,          0,        0;
-a+b*h, 2*a-4*b*h,  -a+3*b*h,          0,        0;
     0,  -a+3*b*h, 2*a-8*b*h,   -a+5*b*h,        0;
     0,         0,  -a+5*b*h, 2*a-12*b*h, -a+7*b*h;
     0,         0,         0,   -a+7*b*h,  a-7*b*h ];

N = @(i,j,x) (x-nodal(j))/(nodal(i)-nodal(j));
ANiipNjjp = @(x, i, ip, j, jp) area(x).*N(i,ip,x).*N(j,jp,x);
M = zeros(n+1,n+1);
for i=1:n
    Vsii =   integral(@(x) ANiipNjjp(x, i, i+1, i, i+1), nodal(i), nodal(i+1));
    Vsiip =  integral(@(x) ANiipNjjp(x, i, i+1, i+1, i), nodal(i), nodal(i+1));
    Vsipip = integral(@(x) ANiipNjjp(x, i+1, i, i+1, i), nodal(i), nodal(i+1));
    M(i:i+1,i:i+1) = M(i:i+1,i:i+1) + rho*[ Vsii, Vsiip; Vsiip, Vsipip ];
end
M

% 解析解との比較
Mp = rho*h/6*[
2*a-b*h,     a-b*h,          0,          0,          0;
  a-b*h, 4*a-8*b*h,    a-3*b*h,          0,          0;
      0,   a-3*b*h, 4*a-16*b*h,    a-5*b*h,          0;
      0,         0,    a-5*b*h, 4*a-24*b*h,    a-7*b*h;
      0,         0,          0,    a-7*b*h, 2*a-15*b*h ];
