Is there such a thing as a matrix of functions? Do we ever put functions as entries of a matrix? If so, are these matrices used in linear algebra or do they have some other special use?
There have been minor not neccessarily conflicts per se, but disagreements on the nature of this question and so I am adding a little statement below to clear this up.
I have noticed that "function" has been interpreted two ways within the answers.


*

*An actual raw function such as merely writing "f". Such a concept is beyond my current understanding (unless I am being stupid somehow), but it is interesting nonetheless.

*A function call returning a value. This is primarily what I meant in my post.
Either one of these is valid. In fact, I think the broadness of this question dictates the fact that people will interpret it differently. In essence, your mileage will vary, and both are so similar from my standpoint that they are all good answers (or good examples if not standalone answers).
 A: A very common example is the matrix of a plane rotation with angle$\theta$ around the origin:
$$\begin{bmatrix}\cos \theta & -\sin \theta\\ \sin \theta & \cos \theta \end{bmatrix}$$                                                      
A: They exist, for example you could have this linear transformation:
$L: \mathbb{P}_2\rightarrow M_{2\times 2}, f\rightarrow L(f):=\begin{pmatrix}
  f'(0) & f(1)  \\
  \int_{-1}^1f(s)\,ds& 0  \\
 \end{pmatrix}$
Functional analysis studies this kind of matrices, I think.
Here $\mathbb{P}_2$ represents the set of all polynomial functions from $\mathbb{R}$ to $\mathbb{R}$, of degree at most 2. The matrix representing this linear transformation is a matrix made of functions appied to some values.
Another example could be the Jacobian matrix, or the Gradient Vector (a vector, but can be seen as a matrix as well)
A: There is the Maurer-Cartan form which is $\mathfrak g$-valued. An example is if we consider the  Lie subgroup $G = SO(2) \subset GL(2,\mathbb R) $, we may parametrize $SO(2)$ by $$g(\theta)= \begin{pmatrix}\cos \theta & -\sin\theta \\\sin \theta & \cos \theta\end{pmatrix}\,\, , \theta \in \mathbb R$$
then the matrix of forms (which are functions) is given by 
$$\omega_G = g^{-1}dg = \begin{pmatrix}0 & -d\theta \\ d \theta & 0\end{pmatrix}$$ 
A: You can define a matrix with elements in any commutative ring, since the only requirement is to be able to perform addition and multiplication with the usual properties.
You even may consider the following $2\times 2$ matrices, with elements that do not belong to the same sets. Such matrices describe the endomorphisms of the direct sum $\;E=U\oplus V$ of two vector spaces $U$ and $V$
$$M=\begin{bmatrix}
f_1&f_2\\g_1&g_2\end{bmatrix},\quad\text{where}\quad\begin{array}{|ll}
f_1\in \mathcal L(U,U),& f_2\in \mathcal L(U,V),\\ g_1\in \mathcal L(V,U),& g_2\in \mathcal L(V,V).
\end{array}$$
You can check one can multiply two such matrices, multiplication of elements being composition of linear maps.
A: One common use of functions in a matrix is the Hessian matrix in multivariable calculus.  This is a matrix of second derivatives with respect to $x_1, x_2, \ldots$.
$$ M = \pmatrix{ \frac{\partial^2 f}{\partial^2 x_1} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial^2 x_2} & \cdots \\ \vdots & \vdots & \ddots
}$$
The eigenvalues of this matrix says a lot about the nature of the function (minimum, maximum, saddle point) and other properties such as stability.
A: Even just the regular first derivative of a function
$$ f:\mathbb{R}^n\rightarrow \mathbb{R}^m$$
Is a matrix of functions, if you define $f$ by its $m$ component scalar functions from $\mathbb{R}^n\rightarrow \mathbb{R}$, or as
$$f(\vec{x})=(f_1(\vec{x}),....,f_m(\vec{x})) $$ with $$
\vec{x}=(x_1,....,x_n)$$ then 
$$
Df(\vec{x})=\begin{bmatrix}\frac{\partial f_1}{\partial x_1}&......&\frac{\partial f_1}{\partial x_n}\\
\frac{\partial f_2}{\partial x_1}&......&\frac{\partial f_2}{\partial x_n}\\
\vdots&&\vdots\\
\frac{\partial f_m}{\partial x_1}&......&\frac{\partial f_m}{\partial x_n}
\end{bmatrix}
$$
Each of the entires above, the partials of the component functions of $f$, are themselves bonafide functions. Or, if you've done some multivariable calculus, the rows are the gradients of the component functions.  
This is what makes linear algebra crucial to studying multivariable calculus. If you want to talk about linear approximations of vector valued functions in something more than the abstract, you have to be comfortable with matrices.  
A: 
Do we ever put functions as entries of a matrix? 

Yes. There have been quite a few fantastic examples given, but I'm not sure they get to the heart of your question.

If so, are these matrices used in linear algebra or do they have some other special use?

The key is not that matrices with functions (or functionals or operators or vectors or matrices etc...) as entries are used in Linear Algebra, it's that we may use Linear Algebra whenever the collection of mathematical objects satisfy certain requirements. 
In linear algebra class, they use numbers as entries in the matrices and vectors, but endeavor to view those matrices and vectors as abstract objects. Using numbers as entries provides context and builds intuition, much like learning arithmetic before algebra.
As for the special uses, I can't think of any common thread between all of the different uses for matrices and vectors with functions as entries. There are quite a few.
A: Here's an example from the study of differential equations:
The Wronskian is the determinant of a matrix in which the $i$-th row represents the $i$-th derivative of $n$ functions. If the determinant is non-zero during a given interval, it demonstrates that the functions are linearly independent along that interval. This is a useful property to know for numerous reasons. 
For example, say you have a second order, linear, homogeneous differential equation, and furthermore say you have found two solutions, $y_1$ and $y_2$, that satisfy the equation. Now let's say you want to find a specific solution given initial conditions. If you can demonstrate that $y_1$ and $y_2$ are linearly independent, then you also know that there exist $c_1$ and $c_2$ such that $y(t) = c_1(y_1) + c_2(y_2)$, given initial conditions $y(t_0) = y_0$ and $y'(t_0) = y_0'$. This is because, when you solve for $c_1$ and $c_2$, the Wronskian ends up in the denominator. If $y_1$ and $y_2$ are linearly dependent (getting back to your question about using linear algebra on matrices with functions), your Wronskian will be zero, and $c_1$ or $c_2$ will not be defined. Thus, linear independence of solutions to a second order, linear, homogeneous differential equation helps demonstrate that a specific solution can be found for a given set of initial conditions.
See here.
