[disp,react]=solve ...

function [disp,react]=solve ...
  (stiff,idu,u,idv,v,idfx,fx,idfy,fy)
%
% [disp,react]=solve(stiff,idu,u,idv, ...
%                    v,idfx,fx,idfy,fy)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function solves for the nodal 
% displacements and reactions.
% 
% stiff   - the stiffness matrix
% idu,u   - nodal indices and magnitudes of all 
%           known horizontal displacements
% idv,v   - nodal indices and magnitudes of all 
%           known vertical displacements
% idfx,fx - nodal indices and magnitudes of all 
%           nonzero known horizontal nodal 
%           forces
% idfy,fy - nodal indices and magintudes of all 
%           nonzero known vertical nodal forces
%
% disp     - computed nodal displacements
% react    - computed nodal forces
%
% User m functions called:  none
%----------------------------------------------

%...Form and solve the system:  STIFF*DISP=F
ns=size(stiff,1); f=zeros(ns,1);
ndc=size([u(:);v(:)]); 
irow=[2*idu(:)-1;2*idv(:)];
%...Form vector f involving known forces and 
%...displacements
f([2*idfx(:)-1;2*idfy(:)])=[fx(:);fy(:)]; 
f(irow)=[u(:);v(:)];
%...Modify rows of the stiffness matrix 
%...corresponding to known displacements
diagv=zeros(ns,1); diagv(irow)=ones(ndc,1); 
s=stiff; s(irow,:)=zeros(ndc,ns); 
s=s+diag(diagv);
%...Solve for displacements and nodal forces
disp=s\f; react=stiff*disp;