how can i write a matlab code for beam vibration and get the natural frequencies and mode shapes

i have a reinforced concrete beam that would be discretized into small elements to account for cracks, the equation of motion for the beam vibration for each element is formulated (the famous equation of motion without damping),each equation needs boundary conditions that are expressed as the continuity conditions between the elements, we need the boundary conditions to get the constants of integration.

 note that the solution is analytical not using the FEM method