data <- matrix(scan(),ncol=2,byrow=TRUE)
1790      3929214
1800      5308483
1810      7239881
1820      9638453
1830     12860702
1840     17063353
1850     23191876
1860     31443321
1870     38558371
1880     50189209
1890     62979766
1900     76212168
1910     92228496
1920    106021537
1930    123202624
1940    132164569
1950    151325798
1960    179323175
1970    203302031
1980    226542199
1990    248709873
2000    281421906

data
nonlin1 <- data.frame(data)
colnames(nonlin1) <- c("year","population")

rm(data)

nonlin1$popmill <- nonlin1$population/1000000 
nonlin1$tyear   <- (nonlin1$year -1790)/10 

# create a function to contain the model and derivatives
popfct  <-  function(b1,b2,b3,x)
  {
     model.func <- b1/(1 +exp(b2 +b3*x))
     derb1      <-  1/(1 +exp(b2 +b3*x))
     derb2      <-  -b1*exp(b2 +b3*x)/(1 +exp(b2 +b3*x))^2
     derb3      <-  -b1*x*exp(b2 +b3*x)/(1 +exp(b2 +b3*x))^2    
     grad.func   <-  cbind(derb1,derb2,derb3)
     attr(model.func, "gradient") <- grad.func
     model.func
  }

popnls <- nls(popmill ~ popfct(b1,b2,b3,tyear), data=nonlin1,
            start=list(b1=350,b2=4.5,b3=-0.3), trace=T  )

summary(popnls)