function direct (ndfree, a, b, c, drp, r, dr, d2r, p, ...
                       omega, dt, nsteps, iprint, nbc, ibc)
%     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%           step by step integration of matrix equations:
%             a*d2r(t)/dt2 + b*dr(t)/dt + c*r(t) = p(t)
%     from initial conditions: r(0), dr(0)/dt, to linear system 
%         s*r(t) = f(t) at each time. (by the wilson method)
%     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%     allocate   a(ndfree,ndfree), b(ndfree,ndfree), ibc(nbc), ...
%                r(ndfree), dr(ndfree), d2r(ndfree), p(ndfree),...
%                drp(ndfree), c(ndfree,ndfree)
      zero = 0. ;
%     initial values of r and dr are passed thru arguments
%     a        = mass (or capacitor) matrix
%     b        = damping (or resistor) matrix
%     c        = stiffness (or inductor) matrix
%     drp      = value of dr at time = t + delt = tau
%     dt       = time step size
%     f        = linear system net forcing, stored over p
%     ibc      = array containing the nbc dof's with zero bc
%     iprint   = no. of integration steps between printing
%     nbc      = no. dof with specified values of zero (if any)
%     nsteps   = no. of integration steps, t_final = dt*nsteps
%     omega    = 1.25 is suggested extrapolation amount
%     p        = forcing function
%     r,dr,d2r = 0,1,2 order deriv. of r wrt t at time=t
%     s        = linear system square matrix, stored over b
%     tau      = extrapolated time step, omega*dt
%               ** initial  calculations **
      disp ('direct step by step integration')
      if  iprint > 1 
        iprint = 1 ;        % data check
      end
      nsteps = (nsteps/iprint)*iprint ;
      if  nsteps <= 0 
        nsteps = iprint ;  % data check
      end % if validation
      delt  = (omega - 1.)*dt  ;  % extrapolation increment
      tau   = omega*dt         ;  % extrapolation amount   
      if  nbc < 1  
        disp ('warning: no constraints in direct \n')
      end % if warning
%      initialize arrays for loads< velocity< acceleration
      p = zeros(ndfree,1); drp = zeros(ndfree,1); d2r = zeros(ndfree,1);
%      overwrite "b" with "s" matrix to save storage
      b = b + 2.*a/tau + tau*c/3. ;
%             ** apply boundary conditions ( to s ) **
      if  nbc > 0 % then apply essential boundary conditions
        for i = 1:nbc
          [b, p] = modful (ndfree, ibc(i), zero, b, p)
        end % do on essential boundary conditions
      end % if boundary conditions needed
%                  ** invert b (once) **
      b = inv (b);  % actually, factor if s is large
%          (approximate the initial value of d2r)
%            ** calculate solution at time t **
      icount = iprint - 2 ;
      for istep = 0: nsteps   % step zero is for initial acceleration
        icount = icount + 1 ;
        t      = dt*(istep-1) ;
        tplus  = t + delt ;              % extrapolate for the load
        if  istep <= 1 
          tplus = zero ;                 % unless initial condition
        end % if
        if  icount == iprint 
          fprintf(1,'\n step number = %5i  time = %e \n', istep, t);
          fprintf(1,'\n   i      r(i)         dr/dt         d2r/dt2 \n')
        end % if
%        forcer is a subr to define the forcing function p(t)
        p = forcer (tplus, ndfree) ;
%        form modified forcing function, f, at t + delt
        for i = 1:ndfree
          p(i) = p(i) + sum ( a(i,:)*( 2.*dr(:)/tau + d2r(:)) ) ...
          - sum ( c(i,:)*(r(:) + 2.*tau*dr(:)/3. + tau*tau*d2r(:)/6.) );
        end % do  for net forcing vector
%           ** apply boundary conditions ( to f ) **
        if  nbc > 0 
          p( ibc(1:nbc) ) = 0 ;
        end % if bc needed
%        solve for drp at time t+delt, s_inverse*f    
        drp = b*p ;  % forward and back substiture if s is large
%        using data at t-dt and t+delt calculate values at t
        for i = 1:ndfree
          drdt   = (1.-1./omega)*dr(i) + drp(i)/omega  ;
          d2rdt2 = d2r(i)*(1.-2./omega)+2.*(drdt-dr(i))/omega/dt  ;
%          approximate the initial value of d2r, if istep <= 1
          if  istep > 1 
            r(i)  = r(i) + dt*(2.*dr(i)+drdt)/3. + dt*dt*d2r(i)/6. ;
            dr(i) = drdt ;
          end % if it is a starter step
          d2r(i) = d2rdt2 ;
        end % do for each degree of freedom
%        output results for time t
        if  icount == iprint % then print this step
          for k = 1:ndfree
            fprintf (1,'%5i %e %e %e \n', k, r(k), dr(k), d2r(k) )
          end % for dof result listing at this time
          icount = 0 ;
        end % if print increment
      end % do over all times
      end % function direct