Inconsistent results of fit_deriv() for different (kind of) models

### Description

While writing a wrapper, to use Astropy (AP) FittableModels in BayesicFitting (BF) I discovered that the method fit_deriv() gives a different result when used on a Gaussian1D and on a Hermite1D.

See python section for some code. According to the code, gp = gauss1d.fit_deriv() results in a list of 3 ndarrays of length 6, while hp = hermite1d.fit_deriv() yields one ndarray with a shape of (6,3). When one would convert gp to a ndarray: gp = numpy.asarray( gp ) it would have a shape of (3,6).

### Expected behavior

The two formats should be the same. I don't care too much what, as long as it is one of the two already present. Dont choose something altogether different. In AstropyModel I check whether it is a list, then turn that into a ndarray which is in the shape I need it. Otherwise, when not a list, I transpose it.

This was probably never noticed because polynomials are linear in their parameters. They never need partials to be solved.

By the way, the documentation in gauss1d claims the fit_deriv() returns a Jacobian matrix, while the doc of Hermite1d (and all other polynomials???) claim it returns a Vandermonde matrix. Actually, it should return the first partial derivatives to all parameters. This is indeed also called the Jacobian matrix and this Jacobian matrix is the same as the Vandermonde matrix in case of polynomials. But only for ordinary polynomials, For Hermite, Chebyshev etc. it is something different. Still the Jacobian and the first derivatives to the parameters.

I think there is also a documentation issue.

### Actual behavior

### Steps to Reproduce

1. [First Step]
2. [Second Step]
3. [and so on...]

```python
# Put your Python code snippet here.

>>> p = numpy.asarray( [1.2, -0.1, 0.3], dtype=float )
>>> x = numpy.linspace( -2, 2, 6, dtype=float )
>>> gm = modeling.models.Gaussian1D()
>>> gp = gm.fit_deriv( x, *tuple( p ) )
>>> print( gp.__class__, gp[0].shape )
(6,)
>>> print( fmt( gp ) )
[[ 0.000 0.001 0.607 0.249 0.000 0.000]
[ -0.000 -0.018 -2.426 1.662 0.001 0.000]
[ 0.000 0.065 2.426 2.771 0.006 0.000]]
>>> hm = modeling.models.Hermite1D( 2 )
>>> hp = hm.fit_deriv( x, *tuple( p ) )
>>> print( hp.__class__, hp[0].shape )
(3,)
>>> print( fmt( hp ) )
[[ 1.000 -4.000 14.000]
[ 1.000 -2.400 3.760]
[ 1.000 -0.800 -1.360]
[ 1.000 0.800 -1.360]
[ 1.000 2.400 3.760]
[ 1.000 4.000 14.000]]

### System Details

import numpy
import astropy.modeling
from BayesicFitting import formatter as fmt ## a nice way to display matrices

Start with the Gaussian1D and Hermite1D fit_deriv() implementations and the documentation describing their return values. Run the supplied Python example, then compare the derivative shapes and terminology across these models; done means the formats and documentation consistently describe the Jacobian, with regression coverage for the reported case.

Written by the indexing model from the issue text.