function [k]=Beam3Noded_1Dk(X,b,h, E)
%clear all;                                   



syms x z real;

L = X(3)-X(1);

PIM=[1        0       0          0         0          0     ;  ...
     0        1       0          0         0          0     ;  ...
     1       L/2   (L/2)^2     (L/2)^3  (L/2)^4     (L/2)^5 ;  ...
     0        1    2*(L/2)   3*(L/2)^2  4*(L/2)^3  5*(L/2)^4; ...
     1        L       L^2       L^3        L^4        L^5   ;  ...
     0        1      2*L       3*L^2    4*L^3        5*L^4  ];
 
 

Ni = [1 x x^2 x^3 x^4  x^5]*inv(PIM);  % 1x6  * 6x6


%ex = -z  d2w / dx2

%B= -z*[diff(Ni,x,2)]

B= -z*[0 0 2  6*x  12*x^2  20*x^3]*inv(PIM);

k= int( int(B'*E*B*b,x,0,L),z,-h/2,h/2);



   