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


D=E;

syms x z real;
L=X(end)-X(1)

PIM=[1 0 0 0; 0 1 0 0; 1 L L^2 L^3;0 1 2*L 3*L^2];

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

%ex = -z  d2w / dx2

B= -z*diff(Ni,x,2);

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




   