mindist(inibod,inid,dt2,init,inidir,iniv,tol,miter,nplot)

function [mdist,flag] = mindist(inibod,inid,dt2,init,inidir,iniv,tol,miter,nplot)

dbstop if error

global deltat g mass pos2 vel2 nbod ndime nsteps color markersize ...
       xmin xmax ymin ymax

%Computing newmark constants
newm.beta  = 0.25;
newm.gamma = 0.5;
newm.b(1) = 1/(newm.beta*dt2*dt2);
newm.b(2) = -1/(newm.beta*dt2);
newm.b(3) = (1-1/(2*newm.beta));
newm.b(4) = newm.gamma/(newm.beta*dt2);
newm.b(5) = 1-newm.gamma/newm.beta;
newm.b(6) = dt2*(1-newm.gamma/(2*newm.beta));

mdist(1:nbod) = Inf;
if init < deltat; init = deltat; end
inidir = inidir/norm(inidir);
position = recpos(init,deltat,pos2);
velocity = recpos(init,deltat,vel2);
nsteps2 = ceil((nsteps - ceil(init/deltat))*deltat/dt2)-1;

rvelocity = iniv*inidir + velocity(:,inibod);
rposition = inid*inidir + position(:,inibod);
raccelera = zeros(2,1);

iposition = rposition;
iposold = iposition;
vold = rvelocity;
aold = raccelera;
matrix = zeros(ndime,ndime);
force = zeros(ndime,1);

length = max(abs(xmin-xmax),abs(ymin-ymax));
flag = 0;

for istep = 1:nsteps2
    error = 1;
    iiter = 0;
    position = recpos(init+istep*dt2,deltat,pos2);
    while (error > tol && iiter < miter)
        iiter = iiter+1;
        force(:) = 0.0;
        matrix(:,:) = 0.0;
        rhs_vec = -(iposition-rposition)*newm.b(1)-rvelocity*newm.b(2)-raccelera*newm.b(3);
            for idime = 1:ndime
                matrix(idime,idime) = newm.b(1);
                force(idime) = rhs_vec(idime);
            end
            for jbod = 1:nbod
                    vec = position(:,jbod)-iposition(:);
                    dist2 = 0;
                    for idime = 1:ndime
                        dist2 = dist2 + vec(idime)*vec(idime);
                    end
                    for idime = 1:ndime
                        matrix(idime,idime) = matrix(idime,idime) ...
                            + g*mass(jbod)   *dist2^(-3/2);
                        for jdime = 1:ndime
                            matrix(idime,jdime) = matrix(idime,jdime) ...
                                - g*mass(jbod)   *(-3)*vec(idime)*vec(jdime)*dist2^(-5/2);
                        end
                        force(idime) = force(idime) + vec(idime)*g*mass(jbod)/(dist2^1.5);
                    end
            end
        u = matrix\force;
        iposition = iposition+u;
        if (iiter == 1); inirhs = force; end
        error = max(norm(force)/norm(inirhs),norm(iposition-iposold)/norm(iposition-rposition));
        iposold = iposition;
    end
    rvelocity = newm.b(4)*(iposition-rposition) + newm.b(5)*vold + newm.b(6)*aold;
    raccelera = newm.b(1)*(iposition-rposition) + newm.b(2)*vold + newm.b(3)*aold;
    rposition = iposition;
    vold = rvelocity;
    aold = raccelera;
    
    for ibod = 1:nbod
        vec = position(:,ibod)-iposition(:);
        dist2 = 0;
        for idime = 1:ndime
            dist2 = dist2 + vec(idime)*vec(idime);
        end
        mdist(ibod) = min(mdist(ibod),dist2);
        if mdist(ibod) < 1e-6*length
            flag = ibod;
            return
        end
    end
    
    if (rem(istep,nplot) == 0)
        clf
        axis([xmin, xmax, ymin, ymax])
        axis off
        %axis equal
        hold on
        for ibod = 1:nbod
            plot(position(1,ibod),position(2,ibod),'-mo',...
                'LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor',color(:,ibod),...
                'MarkerSize',markersize(ibod));
            hold on
        end
        rcol = [0 0 0];
        rsize = 4;
        plot(rposition(1),rposition(2),'-mo',...
            'LineWidth',2,...
            'MarkerEdgeColor','k',...
            'MarkerFaceColor',rcol,...
            'MarkerSize',rsize);
        hold on
        pause(0.001)
    end
end