Using stiffness method to solve 2d Frame

I want to solve a 2D frame like the question in page 28. I follow all the code but i cant run it. The following are my input data and Frame Element Code. What problem do I have and how to solve it? thank you

Input: Nodes: (x, y)

0 0

0 4

3 5

6 4

6 0

Elements: (Node1 Node2), E, A, I,

1 2 200 50e3 22e-6

2 3 200 50e3 22e-6

3 4 200 50e3 22e-6

4 5 200 50e3 22e-6

BCs (Node_Number dof specified_disp)

1 1 0

1 2 0

5 1 0

5 2 0

Nodal loads (NodeNumber dof specified_load)

2 1 -15

2 2 5

2 3 -25/3

3 1 -30

4 1 -15

4 2 -5

4 3 25/3

Frame element code function [Kel] = FrameElement2D(nodexy, E, A, I)

E1 = [ (nodexy(2,1)-nodexy(1,1))

       (nodexy(2,2)-nodexy(1,2))];

L = norm(E1);

E1 = E1/L;

E2 = [-E1(2) E1(1)];

Kel_bend = [ 12*E*I/(L^3) 6*E*I/(L^2) -12*E*I/(L^3) 6*E*I/(L^2) ;

             6*E*I/(L^2)   4*E*I/L      -6*E*I/(L^2)  2*E*I/L      ;

             -12*E*I/(L/3)  -6*E*I/(L^2) 12*E*I/(L^3)  -6*E*I/(L^2);

             6*E*I/(L^2)    2*E*I/L      -6*E*I/(L^2)  4*E*I/L    ];

Kel_axial = E*A/L*(1 -1; -1 1);

Kel_Loc([1,4],[1,4]) = Kel_axial;

Kel_Loc([2,3,5,6],[2,3,5,6]) = Kel_bend;

Qrot = [E1;E2]; %transforms global to element d_E=Q d_G

Qrot(3,3) = 1;

Tmatrix = [Qrot zeros(3); zeros(3) Qrot];

Kel = Tmatrix'*Kel_Loc*Tmatrix;

Answer 1

[Kevin Díaz]

I've been tryin' to run this same code and still having trouble with it. Have you found any solution?. If you so, please share it, I'd thank you a LOT.