Resampling scattered data with MATLAB's $\tt interp$ and $\tt resample$
In your question, you state that your ultimate goal is to compute the discrete Fourier transform (DFT) from nonuniformly-sampled input data. It is not necessary to resample your data at all.  In fact, there are several fast algorithms, known variously as the unequally spaced fast Fourier transform (USFFT) and the nonuniform fast Fourier transform (NUFFT) that do exactly what you want: take nonuniform samples as input and produce equally spaced samples of the DFT as output. The remainder of your question, dealing with issues related to MATLAB's resampling routines, are a red herring.

Here are some publications describing the USFFT algorithms:


G. Beylkin, On the Fast Fourier Transform of Functions With
Singularities,  Applied and Computational Harmonic Analysis, 2, pp.
363-381, 1995. Link to preprint.
A. Dutt and V. Rokhlin, Fast Fourier Transforms for Nonequispaced Data, SIAM J. Sci. Comput., 14, pp. 1368, 1993.


There are many commercial and free implementations of these algorithms. You can access some of them here.

EDIT: Here is a "quick and dirty" approach that may help you.

Based on your problem description and comments, you are just trying to get a sense of the spectral content of your signal, which does not require high accuracy. Here is one way to formulate your problem so that the spectrum can be easily estimated using any of the USFFT implementations in the above link.

You have a signal $f(t)$ defined on the $T$-duration interval $t\in[-T/2,T/2]$, and you wish to estimate its Fourier transform
$$\hat{f}(\omega)=\int_{-T/2}^{T/2}f(t)e^{-i 2 \pi \omega t}\,dt$$
at a number of values of $\omega$. You have collected $N$ samples of your function
$$f_n=f(t_n),\quad n=1,\ldots,N$$
at the unequally spaced time locations $t_n$, $n=1,\ldots,N$.

In general, you would solve this by formulating a least squares problem for the Fourier transform $\hat{f}$, and using the USFFT to rapidly solve the system of equations, but since you do not require high accuracy, here is a "quick and dirty" solution.

Use the trapezoid rule from Calculus 101 to discretize the above integral, giving the approximation
$$\hat{f}(\omega) \approx \frac{\Delta_1}{2}f_1 e^{-i2\pi\omega t_1}+\left(\sum_{n=2}^{N-1}\frac{\Delta_{n-1}+\Delta_{n}}{2}f_n e^{-i2\pi\omega t_n}\right)+\frac{\Delta_{N-1}}{2}f_N e^{-i2\pi\omega t_N}$$ where
$$\Delta_{n}=t_{n+1}-t_n$$ 
is the temporal spacing between samples $n+1$ and $n$.

Define the new sequence 
$$g_n=\begin{cases}\frac{\Delta_1}{2}f_1, &amp; n=1 \\ \frac{\Delta_{n-1}+\Delta_{n}}{2}f_n, &amp; n=2,\ldots,N-1 \\ \frac{\Delta_{N-1}}{2}f_N, &amp; n=N\end{cases}$$
Now pick some frequency grid you care about. For example, say you want $2000$ equally spaced frequencies between $0$ Hz and $25$ Hz, and define
$$\omega_m = \frac{m-1}{2000-1}25,\quad m=1,\ldots,2000.$$

Now your problem amounts to evaluating the sum
$$\hat{f}(\omega_m) = \sum_{n=1}^N g_n e^{-i2\pi\omega_m t_n},\quad m=1,\ldots,M.$$

This looks like an FFT, except the $t$-values are unequally spaced and the $\omega$-values are unusual. If you do not care about speed, then just evaluate the sum directly using your favorite programming language.  It will take a minute or two of computer time, but so what? If you are just doing this a few times to get a sense of what your signal's spectrum looks like, then it does not matter how fast it is. If at some point in the future you decide that you do care about speed, then the sum can be evaluated extremely rapidly using a single call to a USFFT library, which you can download from the link I posted (or you can use one of the many commercial implementations).  Making the switch to the USFFT will only require you to modify a handful of lines of code.