how to solve this motion equation?

i need to write a function that describes the equations of motin (added in png)

Answer 1

[Torsten]

function main
  y0 = [0 pi/2 0 pi/2];  %[y1 y2 y3 y4]=[theta1 theta1dot theta2 theta2dot]
  tspan = [0 144*pi];
  [T,Y] = ode15s(@fun,tspan,y0);
  plot(T,[Y(:,1),Y(:,3)])
end  

function f = fun(t,y)
  I1 = 1;
  I2 = 1;
  l1 = 1;
  l2 = 1;
  m1 = 2;
  m2 = 3;
  g = 9.8;
  b1 = 0.5;
  b2 = 0.2;
  A = zeros(2,2);
  b = zeros(2,1);
  A(1,1) = I1+(m1+4*m2)*l1^2;
  A(1,2) = 2*m2*l1*l2*cos(y(1)-y(3));
  A(2,1) = 2*m2*l1*l2*cos(y(1)-y(3));
  A(2,2) = I2+m2*l2^2;
  b(1) = -(2*m2*l1*l2*sin(y(1)-y(3))*y(4)^2 + ...
         (m1+2*m2)*l1*g*sin(y(1)) + ...
         b2*(y(2)-y(4)) + b1*y(2));
  b(2) = -(-2*m2*l1*l2*sin(y(1)-y(3))*y(2)^2 + ...
         m2*l2*g*sin(y(3)) - b2*(y(2)-y(4)));
  sol = A\b;
  f = [y(2);sol(1);y(4);sol(2)];
end