function X = stfew(t,x)

global GK GM evec L B ic voltage v_time FM Flag freq ...
vdc vac omega gap e_force_elems ncnt

ic = ic +  1;
v_time(ic) = t;
volt = vdc + vac*sin(omega*t);
voltage(ic) = volt;

N = size(GK,2);

% Collect only modal coordinates skipping their derivatives
idvec = 1:2:2*FM;
q = x(idvec);

% Get back original coordinates (nodal disps & thetas)
posnow = evec(:,1:FM)*q;

% Collect only displacements skipping thetas
idvec = 1:2:size(GK,2);
ynow = posnow(idvec);

% Since no forces exist on theta state variables, zeros
% are included in between to make up the full force vector
forcevec = zeros(2*size(ynow,1),1);

j = 0;
for i = 1:size(e_force_elems,2),
    elemID = e_force_elems(i);
    node1ID = ncnt(elemID,1)-1; dof1 = 2*(node1ID-1) + 1;
    node2ID = ncnt(elemID,2)-1; dof2 = 2*(node2ID-1) + 1;
% In computing the dof number, the boundary condition
% is assumed to be fixed-fixed.
%   To convert force to micron  units,  mulitply by 1E6
    mag = 8.854e-12 * L(elemID)*B(elemID) * volt^2 *1e6/ 2;
%   Fringing field correction
    yave = ( ynow(node1ID) + ynow(node2ID) ) / 2;
    ffc = (1+0.65* (gap-yave)/B(elemID));
    forcevec(dof1) = forcevec(dof1) + ffc*mag / (gap - yave)^2 / 2;
    forcevec(dof2) = forcevec(dof2) + ffc*mag / (gap - yave)^2 / 2;
end

% Convert forcevec to modal space
eigforce = evec'*forcevec;

% Compute derivatives to be supplied to ODE23
j = 0;
for i = 1:FM,
  disp = x(j+1);
  vel = x(j+2);
  j = j + 1; X(j) = x(j+1);
  j = j + 1; X(j) = -GK(i,i)/GM(i,i)*disp + eigforce(i)/GM(i,i) ;
%...                    - damp*vel/GM(i,i);
end

X = X';

%clc home
%t
%volt
%mag
%eigforce
%maxdisp = max(abs(ynow))
%maxeig = max(eigforce)








