four-bar mechanism full analysis (Double-Rocker)

Hello,

I am a self-learner beginner in MATLAB and I am working on a full analysis of a 4-bar mechanism for a project at my university .

I started with a crank rocker that I wanted to turn into a Double rocker but by doing so the code doesn't work .

I was trying to do those three things and I am asking for your help :

1)I wanted to make a plot between theta input (2) and theta output (4) but i wanted to break the loop so th2 starts oscilating from 30* to 150* ccw then cw .

2)I was making a force analysis but I don't know how to plot max power.

The values of angular velocities and acceleration are just for demonstration , they will be changed to their real values after the syntax is done.

------------------

Here is the syntax :

clc; clear all;

t=0:0.05:10;
 
l1=5; %fixed
l2=2; %input link
l3=3; %coupler
l4=4; %driver
w2 = 1 ; %ccw 
a2= 1; %ccw
m= 2;
g=9.81;


% Plotting

for th2=0:0.1:2*pi
    
   %position 
   O= [0;0];
   Q= l1*[1;0];
   A=[l2*cos(th2); l2*sin(th2)];
   lD=sqrt(l2^2 + l1^2 - 2*l2*l1*cos(th2));
   alpha=asin(l2*sin(th2)./lD);
   beta=acos((lD.^2+l4^2-l3^2)./(2*lD*l4));
   B=[l1-l4*cos(alpha+beta);l4*sin(alpha+beta)];
   th3= atan((B(2) - A(2))/(B(1) - A(1)));
   th4= atan(B(2)/(B(1)-l1));
   
   
   pointO=viscircles(O,0.01);
   pointA=viscircles(A,0.01);
   pointB=viscircles(B,0.01);
   pointQ=viscircles(Q,0.01);
   
   
   l2_bar = line([O(1) A(1)],[O(2) A(2)]);
   l3_bar = line([l2*cos(th2) l2*cos(th2)+l3*cos(th3)],[l2*sin(th2) l2*sin(th2)+l3*sin(th3)],'linewidth',0.5);
   l4_bar = line([l1 l1+l4*cos(th4)],[0 l4*sin(th4)],'linewidth',2);
   
   
   
   %velocity
   VA=[-w2*l2*sin(th2) w2*l2*cos(th2)];
   w3=([VA(1)*cos(th4)+VA(2)*sin(th4)]/l3*sin(th3-th4));
   w3=([-VA(1)*cos(th3)-VA(2)*sin(th3)]/l3*sin(th4-th3));
   VBA=[-w3*l3*sin(th3) w3*l3*cos(th3)];
   VB=VA+VBA;
   
   %acceleration
   AA=[l2*(-w2.^2*cos(th2)-a2*sin(th2))  l2*(-w2.^2*sin(th2)+a2*cos(th2))];
   a3=(([(AA(1)*cos(th4))+(AA(2)*sin(th4))-(l3*w3.^2*cos(th3-th4))]/l3*w3.^2*sin(th3-th4)));
   a4=([(-AA(1)*cos(th3))-(AA(2)*sin(th3))-(l4*w4.^2*cos(th4-th3))]/l4*w4.^2*sin(th4-th3));
   ABA=[l3*(-w2.^2*cos(th3)-a2*sin(th3))  l3*(-w3.^2*sin(th)+a2*cos(th3))];
   AB=AA+ABA;
   
   
  
   
   
   
   
    pause(0.01)
    
    delete(pointO);
    delete(pointA);
    delete(pointB);
    delete(pointQ);
    delete(l2_bar);
    delete(l3_bar);
    delete(l4_bar);
   
   axis(gca,'equal');
 
  axis([-3 8 -5 5]);
end