Ok so , divide the beam or beams into sections between nodes where you can solve displacement and rotation of each node . Set them up equal length from each other because it makes it much easier . For each section you have a local stiffness matrix EIL 3 . 12 6 L 12 6 L 6 L 4 L 2 6 L 2 L 2 12 6 L 12 6 L 6 L 2 L 2 6 L 4 L 2 where columns from left to right describe displacement then rotation of left node then displacement then rotation of right node . Create one for each section . I lie to create a multidimensional matrix kzeros 4 , 4 , % number of elements %. If beam remains constant EIL 3 for each section you can multiply that later in a global stiffness matrix where the columns left to right are displacement then rotation of each node from first to last . In the global stiffness matrix I like to call K , add the shared nodes for different elements in the respective locations within the global matrix . the EulerBernoulli theorem for the global matrix goesF 1 y M 1 F 2 yM 2 F 3 yM 3 . V 1 phi 1 V 2 phi 2 V 3 phi 3 . The left side of the equation is the applied loads and moments you can create a matrix for that . Now solving the matricies is where I have trouble not knowing much linear algebra for MATLAB . You can try linsolve , something like FK or vise versa .