function psip=nec4bar(t,psi)
% sistema não linear, excitação (momento) constante, c.i. não nulas

global L0 L1 L2 L3

% mecanismo de 4 barras
  % distâncias dos centros de massa (m)
Lc1=L1/2; Lc2=L2/2; Lc3=L3/2;
  % massas das barras (kg)
m1=1.0; m2=1.0; m3=1.0;
  % momentos de inércia de massa das barras (kg.m2)
I1=m1*L1^2/12; I2=m2*L2^2/12; I3=m3*L3^2/12;

% torque de entrada (N.m)
Tt=6;

% aceleração da gravidade (m/s2)
g=9.81;

% grandezas cinemáticas
  % variável primária - barra 1
teta=psi(1);
  % variáveis secundárias
    % barra 2
k1=-2*L1*L3*sin(teta);
k2=2*L3*(L0-L1*cos(teta));
k3=L0^2+L1^2-L2^2+L3^2-2*L0*L1*cos(teta);
fi=wrapToPi(2*atan2(-k1-sqrt(k1^2+k2^2-k3^2),k3-k2));
    % barra 3
alfa=atan2(-L1*sin(teta)+L3*sin(fi),L0-L1*cos(teta)+L3*cos(fi));

% coeficientes intermediários
J1=0.5*(m1*Lc1^2+I1+m2*L1^2);
J2=0.5*(m2*Lc2^2+I2);
J3=0.5*(m3*Lc3^2+I3);
S1=(L1*sin(fi-teta))/(L2*sin(alfa-fi));
S2=(L1*sin(alfa-teta))/(L3*sin(alfa-fi));
P1=m2*L1*Lc2;
C1=cos(teta-alfa);

% derivadas intermediárias
dS1dte=-(L1*cos(fi-teta))/(L2*sin(alfa-fi));
dS1dal=-(L1*sin(fi-teta)*cos(alfa-fi))/(L2*sin(alfa-fi)^2);
dS1dfi=-(2*L1*sin(alfa-teta))/(-L2+L2*cos(2*alfa-2*fi));
dS2dte=-(L1*cos(alfa-teta))/(L3*sin(alfa-fi));
dS2dal=(2*L1*sin(fi-teta))/(-L3+L3*cos(2*alfa-2*fi));
dS2dfi=(L1*sin(alfa-teta)*cos(alfa-fi))/(L3*sin(alfa-fi)^2);
dC1dte=sin(alfa-teta);
dC1dal=-sin(alfa-teta);
dGdte=-(m1*Lc1+m2*L1)*g*cos(teta);
dGdal=-m2*g*Lc2*cos(alfa);
dGdfi=-m3*g*Lc3*cos(fi);

% coeficientes em forma compacta
M=2*(J1+J2*S1^2+J3*S2^2+P1*C1*S1);
N=(2*J2*S1*(dS1dte+S1*dS1dal+S2*dS1dfi)+...
   2*J3*S2*(dS2dte+S1*dS2dal+S2*dS2dfi)+...
   P1*(C1*(dS1dte+S1*dS1dal+S2*dS1dfi)+S1*(dC1dte+S1*dC1dal)));
R=-dGdte-S1*dGdal-S2*dGdfi;

psip=[psi(2); -N/M*psi(2)-R/M+Tt/M];