MATLAB: Approximate tomorrow's temperature with 2nd, 3rd and 4th polynomial using the Least Squares method. The following is Exercise 3 of a Numerical Analysis task I have to do as part of my university course on the subject.

Find an approximation of tomorrow's temperature based on the last 23
  values of hourly temperature of your city ( Meteorological history for
  Thessaloniki {The city of my univ} can be found here:
  http://freemeteo.com)
You will approximate the temperature function with a polynomial of
  2nd, 3rd and 4th degree, using the Least Squares method. Following
  that, you will find the value of the function at the point that
  interests you. Compare your approximations qualitatively and make a
  note to the time and date you're doing the approximation on.

Maybe it's due to fatigue due to doing the first two tasks without break, or it's my lack of experience on numerical analysis, but I am completely stumped. I do not even know where to start.
I know it's disgusting to ask for a solution without even showing signs of effort, but I would appreciate anything. Leads, tutorials, outlines of the things I need to work on, one after the other, anything.
I'd be very much obliged to you.
NOTE: I am not able to use any MATLAB in-built approximation functions.

Let $\mathbf{t}\in\mathbb{R}^{23}$ be the last 23 samples you have. To fit these to an $N^\mathrm{th}$ order  polynomial in terms of the hour, i.e., $t = \sum_i p_i h^i$, where $\mathbf{p}\in\mathbb{R}^{N}$ is the vector of coefficients and $h\in\mathbb{R}$ is the time in hours, first set up the system using the measurements you have and the times at which they were taken. According to the problem, the times were taken at $h = 0,1,\dots$, so let us define $h_i$ to be equal to the value $i$. Then
$$
\left[
\begin{array}{ccccc}
h_0^0 & h_0^1 & h_0^2 & h_0^3 & \cdots \\
h_1^0 & h_1^1 & h_1^2 & h_1^3 & \cdots \\
\vdots & \vdots & & \ddots
\end{array}
\right]
\mathbf{p} = \mathbf{H} \mathbf{p} = \mathbf{t}
$$
If you write the above out in long form, you'll see that all it says is that
$$
p_0 \cdot h_i^0 + p_1 \cdot h_i^1 + p_2 \cdot h_i^3 + \dots = t_i
$$
Now to solve for $\mathbf{p}$, you just use the normal equation:
$$
\mathbf{p} = \left( \mathbf{H}^{\mathrm{T}} \mathbf{H} \right)^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{t}
$$
In Matlab, setting up and solving for the above is straightforward. Once you have $\mathbf{H}$ and $\mathbf{t}$ defined,
p = (H' * H) \ H' * t;

