function W = ComputeWorkDone(p,i)
W = zeros(length(p.xdot),1);
for ii = i;
  for j = 1:p.dof
    %j = 3;
    W = W + p.dt*cumtrapz(p.Forces{ii}(:,j).*p.xdot(:,j));
  end
end