Here is the corrected code.
Note the errors with:
(1) symbolic math expression and solving the equations: syms, solve()
(2) singularity problem with the initial angle at phi=0 degree (0.01*pi taken to avoid this issue). To speed up simulation, take a larger step for phi, e.g.: phi= =0.01*pi:pi/10:6*pi
(3) Suggestion: instead of displaying all values with fprintf(), try to store them.
Have fun!
clear all; clc; close all
% Input data
AB=0.22; %(m)
AC=0.08; %(m)
CD=0.20; %(m)
DE=0.60; %(m)
M = moviein(12);
incr = 0;
VB23=zeros(1,360);
VB3=zeros(1,360);
VBD=zeros(1,360);
VBE=zeros(1,360);
AB1=zeros(1,360);
AB23=zeros(1,360);
AB3=zeros(1,360);
ABD=zeros(1,360);
ABE=zeros(1,360);
TANPHI=zeros(1,360);
i=1;
for phi =0.01*pi:pi/90:6*pi
%phi = 4*pi/3 ; %(rad)
fprintf('\nphi =%g (rad)\n', phi)
n = 500*2 * pi/60;
omega1 = [0 0 n ]; % (rad/s)
fprintf('omega1 = omega2 = [ %g, %g, %g ] (rad/s)\n', omega1)
% Position of joint A
fprintf('\nPoint A \n')
xA=0;yA=0;rA = [xA yA 0];
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rA = [ %g, %g, %g ] (m)\n', rA)
vA=[0 0 0];%(m/s) % velocity of A (fixed)
fprintf('vA = [ %g, %g, %g ] (m/s)\n', vA)
fprintf('|vA|= %g (m/s)\n',norm(vA))
% Position of joint B1
fprintf('\nPoint B \n')
xB=AB*cos(phi); yB = AB*sin(phi); rB = [xB yB 0];
fprintf('rB = [ %g, %g, %g ] (m)\n', rB)
vB1 = vA + cross(omega1,rB); % velocity of B1
fprintf('vB1 = vB2 = [ %g, %g, %g ] (m/s)\n', vB1)
fprintf('|vB1| = |vB2| = %g (m/s)\n', norm(vB1))
vB2x=vB1(1);
vB2y=vB1(2);
omega11=norm(omega1);
aB1 = [-(omega11)^2*xB -(omega11)^2*yB 0];
fprintf('aB1 = aB2 = [ %g, %g, %g ] (m/s2)\n', aB1)
%%%%%%
% tan phi3
tanphi3 = (yB-yC)/(xB-xC);
phi3=atan(tanphi3);
fprintf('phi 3 = %g (rad/s)\n', phi3)
cosphi3= cos(phi3);
fprintf('cosphi3 = %g (rad/s)\n', cosphi3)
sinphi3= sin(phi3);
fprintf('sinphi3 = %g (rad/s)\n', sinphi3)
%%%
syms vB23x vB3x
eqnB3x = vB2x + vB23x == vB3x ;
eqnB23x = -vB3x*((xC-xB)/(yC-yB))==vB2y + vB23x*((yB-yC)/(xB-xC));
solB = solve(eqnB3x, eqnB23x, vB3x,vB23x);
vB3x= eval(solB.vB3x);
vB23x= eval(solB.vB23x);
syms vB3y vB23y
eqnB3y = vB2y + vB23y == vB3y ;
eqnB23y = -vB3y*((yC-yB)/(xC-xB))==vB2y + vB23y*((xB-xC)/(yB-yC));
solB = solve(eqnB3y, eqnB23y, vB3y,vB23y);
vB3y= eval(solB.vB3y);
vB23y= eval(solB.vB23y);
vB3 = [vB3x vB3y 0];
vB23 = [vB23x vB23y 0];
omega3 = [0 0 vB3y/(xB-xC)];
fprintf('omega3 = [ %g, %g, %g ] (rad/s)\n', omega3)
fprintf('vB3 = [ %g, %g, %g ] \n', vB3)
fprintf('|vB3| = %g (m/s)\n', norm(vB3))
fprintf('vB23 = [ %g, %g, %g ] \n', vB23)
syms alpha3 aB23
eqnalpha3 = -alpha3*(yC-yB) == aB1(1) + aB23*cos(phi3) -vB23(1)*2*sin(90*pi/360-phi3);
eqnaB23 = alpha3*(xC-xB) == aB1(2) + aB23*sin(phi3) +vB23(2)*2*cos(90*pi/360-phi3);
solalpha = solve(eqnalpha3, eqnaB23, alpha3, aB23);%giai pt
alpha3pos = eval(solalpha.alpha3);
aB23pos = eval(solalpha.aB23);
fprintf('alpha3 = %g \n', alpha3pos )
fprintf('aB32pos = %g \n', aB23pos )
aB23p = [aB23pos*sinphi3 aB23pos*cosphi3 0];
fprintf('aB32 = [ %g, %g, %g ] \n', aB23p)
aB3 = alpha3pos*(rB-rC);
fprintf('aB3 = [ %g, %g, %g ] \n', aB3)
% Position of joint C
fprintf('\nPoint C \n')
xC=0.08; yC = 0; rC = [xC yC 0];
fprintf('rC = [ %g, %g, %g ] (m)\n', rC)
vC = [0 0 0];
fprintf('vC = 0 (m/s)\n')
fprintf('|vC|= 0\n')
% Position of joint D
fprintf('\nPoint D \n')
syms xDsol yDsol
eqnD1 = (xDsol - xC)^2 + (yDsol - yC)^2 == CD^2;
eqnD2 = (yB/(xB-xC))*(yDsol/(xDsol - xC)) == -1;
solD = solve(eqnD1, eqnD2, xDsol, yDsol);%giai pt
xDpositions = eval(solD.xDsol);
yDpositions = eval(solD.yDsol);
xD1= xDpositions(1);
xD2= xDpositions(2);
yD1= yDpositions(1);
yD2= yDpositions(2);
xD = xD1;
yD=yD1;
rD = [xD yD 0];
fprintf('rD = [ %g, %g, %g ] (m)\n', rD)
vD = vC + cross(omega3, rD-rC);
fprintf('vD = vD5 = [ %g, %g, %g ] (m/s)\n', vD)
fprintf('|vD| = %g (m/s)\n', norm(vD))
vDx=vD(1);
vDy=vD(2);
aD= alpha3pos*(rD-rC) - norm(omega3)^2*(rD-rC);
fprintf('aD = [ %g, %g, %g ] (m/s2)\n', aD)
% Position of joint E
syms xEsol
eqnE = (xEsol - xD)^2 + (yD)^2 == DE^2 ;
solE = solve(eqnE, xEsol);
xEpositions = eval(solE);
xE1= xEpositions(1);
xE2= xEpositions(2);
if xE1> xC
xE = xE1;
else
xE = xE2;
end
yE=0;
rE = [xE yE 0];
fprintf('\nPoint E \n')
fprintf('rE = [ %g, %g, %g ] (m)\n', rE)
syms vEx omega4
eqnE1 = vDx - omega4*(yE-yD) == vEx ;
eqnE2 = vDy + omega4*(xE-xD) == 0;
solD = solve(eqnE1, eqnE2, vEx, omega4);%giai pt
vEx = eval(solD.vEx);
omega4x = eval(solD.omega4);
vE = [vEx 0 0];
omega4 = [0 0 omega4x];
fprintf('vE = [ %g, %g, %g ] (m/s)\n', vE)
fprintf('|vE| = %g (m/s)\n', norm(vE))
fprintf('omega4 = [ %g, %g, %g ] \n', omega4)
aDx=aD(1);
aDy=aD(2);
syms aEx alpha4
eqnaEx = aEx == aDx + alpha4*(yE-yD) - (norm(omega4))^2*( yE-yD);
eqnalpha4 = aDy+alpha4*(xE-xD) - (norm(omega4))^2*(xE-xD)==0;
solaE = solve(eqnaEx, eqnalpha4, aEx, alpha4);%giai pt
aExpos = eval(solaE.aEx);
alpha4pos = eval(solaE.alpha4);
fprintf('alpha4 = %g \n', alpha4pos )
fprintf('aE = [%g 0 0] \n', aExpos )
% Graphic of the mechanism
figure(1)
plot([xA,xB],[yA,yB],'k-*',...
[xB,xC],[yB,yC],'b-*',...
[xC,xD],[yC,yD],'g-*',...
[xD,xE],[yD,yE],'r-*','LineWidth',1.5)
% adds major grid lines to the current axes
%grid on,...
xlabel('x (m)'), ylabel('y (m)'),...
title('positions for \phi = 120 (deg)'),...
text(xA,yA,'\leftarrow A = ground',...
'HorizontalAlignment','right'),...
text(xB,yB,'--B'),...
text(xC,yC,'--C'),...
text(xD,yD,'--D'),...
text(xE,yE,'--E'),grid;
%axis([-0.3 0.8 -0.3 0.45])
xlim([-0.3 1.3]);
ylim([-0.3 0.9]);
% end of program
%incr = incr + 1;
%M( : , incr) = getframe; % record the movie
VB23(i)=norm(vB23);
VB3(i)=norm(vB3);
VD(i)=norm(vD);
VE(i)=norm(vE);
AB1(i)=norm(aB1);
AB23(i)=norm(aB23p);
AB3(i)=norm(aB3);
AD(i)=norm(aD);
AE(i)=norm(aExpos);
i=i+1;
end % end for
figure(4)
subplot ( 2, 2,1)
plot(VB23,'r');
title('|vB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(VB3,'g');
title('|vB3|')
subplot ( 2, 2,3)
plot(VD,'k');
title('|vD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(VE,'b');
legend('y = sin(x)','y = cos(x)');
title('|vE|')
xlabel('0 < phi < 6*pi') % x-axis label
figure(5)
subplot ( 2, 2,1)
plot(AB23);
title('|aB23|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,2)
plot(AB3);
title('|aB3|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,3)
plot(AD);
title('|aD|')
xlabel('0 < phi < 6*pi') % x-axis label
subplot ( 2, 2,4)
plot(AE);
title('|aE|')
xlabel('0 < phi < 6*pi') % x-axis label
%movie2avi(M,'exam7_9_velocity.avi');
