Use par(mfrow=c(2,3)) to make the following plot be arranged in one 2 * 3 grid.
If you want fine control, keep reading here(layout), here(ggplot+gridExtra) 
png(filename="C:\\Users\\datafireball.com\\Documents\\R\\stackoverflow_7144118.png")
par(mfrow=c(3,2))
plot(rnorm(10), rnorm(10))
plot(rnorm(10), rnorm(10))
plot(rnorm(10), rnorm(10))
plot(rnorm(10), rnorm(10))
plot(rnorm(10), rnorm(10))
plot(rnorm(10), rnorm(10))
dev.off()
You can remove the first and last line so you can print it to the standout output. 

Update: In your case, looks like par(mfrow) won't work, because I don't think it is actually calling the base plot method, instead, the return object from the fitMod method is actually a type called "trellis", which belongs to the lattice package. If you want to know more about trellis, read here. However, if you just want to know how to get it work, I got it working with the grid.arrange method from gridExtra. 
library(DoseFinding)
library(gridExtra)
data(biom)
# here, the bnds argument has been ignored so the default value from defBnds will be applied. 
fitemax <- fitMod(dose, resp, data=biom, model="emax")
p1 <- plot(fitemax)
fitlinearlog <- fitMod(dose, resp, data=biom, model="linlog")
p2 <- plot(fitlinearlog)
fitlinear <- fitMod(dose, resp, data=biom, model="linear")
p3 <- plot(fitlinear)
fitquadratic <- fitMod(dose, resp, data=biom, model="quadratic")
p4 <- plot(fitquadratic)
fitexponential <- fitMod(dose, resp, data=biom, model="exponential")
p5 <- plot(fitexponential)
fitlogistic <- fitMod(dose, resp, data=biom, model="logistic")
p6 <- plot(fitlogistic)
grid.arrange(p1, p2, p3, p4, p5, p6)
# Message: Need bounds in "bnds" for nonlinear models, using default bounds from "defBnds".

Is this the output you want?