To do it more accurately, you need to pick ellipsoid parameters, which measure the distortion of the oblate spheroid. For North America, the Clarke 1866 parameters are often used. The Clarke 1866 parameters are:

      $a = 6378206.4;
      $b = 6356583.8;
[download]

$a corresponds to a North/South radius, and $b corresponds to an East/West radius. Then you can compute the earth radius at a given latitude. I only have the calculation in C++ code at the moment (Latitude and LatLon are classes holding the numeric values):

// Reduced or Parametric Latitude
   Latitude Ellipsoid::
Eta(const Latitude & lat) const {
   // NOTE: sqrt(1-e^2) == b/a
   return Latitude(atan(b/a*tan(lat*deg2rad))*rad2deg);
}  // Eta


// Radius of Ellipsoid at Latitude.
   double Ellipsoid::
Radius(const Latitude & lat) const {
   Latitude eta = Eta(lat);
   double x = a * cos(eta*deg2rad);
   double y = b * sin(eta*deg2rad);
   return sqrt(x*x + y*y);
}  // Radius
[download]

Then you can compute a new location from a lat/lon point and an azimuth:

// Given a Lat,Lon point, a distance (in meters), and an azimuth (in d
+egrees),
// this routine returns a new Lat,Lon point.
   LatLon Ellipsoid::
NewLatLon(const LatLon & lat_lon, double c, double Az) const {
   // Convert to radians
   double lat = lat_lon.lat * deg2rad;
   // We do not need lon in radians.
   Az *= deg2rad;

   // Equations (5-5) and (5-6) in "Map Projections--A Working Manual"
+, Page 31
   double cosAz = cos(Az);
   double sinAz = sin(Az);
   double sinlat = sin(lat);
   double coslat = cos(lat);

   c /= Radius(lat_lon.lat);
   double cosc = cos(c);
   double sinc = sin(c);

   double y = sinc*sinAz;
   double x = coslat*cosc - sinlat*sinc*cosAz;
   double at = (x == 0.0  &&  y == 0.0) ? 0.0 : atan2(y,x);
   lat = asin(sinlat*cosc + coslat*sinc*cosAz) * rad2deg;
   double lon = lat_lon.lon + at * rad2deg;

   return LatLon(lat, lon);
}  // NewLatLon
[download]