Another way (by Euler, I think), from the geometric sum:

$$1 + x + x^2 + cdots + x^n = frac{x^{n+1}-1}{x-1} tag{1}$$

Differentiate both sides and multiply by $x$:

$$x + 2 x^2 + 3 x^3 + cdots + n x^{n} = frac{n x^{n+2}-(n+1) x^{n+1} +x}{(x-1)^2} tag{2}$$

Differentiating once more, we get on the LHS

$$1 + 2^2 x + 3^2 x^2 + cdots + n^2 x^{n-1} tag{3}$$

which, evaluated at $x=1$ gives our desired sum $sum_{k=1}^n k^2$. Hence, we just need to compute the derivative on the RHS in $(2)$ (eg, with L'Hopital rule) and evaluate it at $x to 1$.

It should be evident that this procedure also can be applied (though it also turns more cumbersome) for sums of higher powers.

(Update) here's a concrete computation, by doing a series expansion of around $(x-1)$, applying the binomial theorem on the RHS of $(1)$. Let

$$begin{align}
g(x)&=frac{x^{n+1}-1}{x-1}\
&=frac{left(1+(x-1)right)^{n+1}-1}{x-1}\
&={n+1 choose 1}+{n+1 choose 2}(x-1)+{n+1 choose 3}(x-1)^2+Oleft((x-1)^3right) tag{4}
end{align}$$

Deriving, multiplying by $x$ and deriving again: $$(g'(x) , x)'={n+1 choose 2}+{n+1 choose 3}2 , x + O(x-1) tag{5}$$
which evaluated at $x=1$ gives the desired answer:
$${n+1 choose 2}+2{n+1 choose 3} =frac{n(n+1)(2n+1)}{6} $$

We can apply the same procedure for higher powers. For example:
$$ sum_{k=1}^n k^3={n+1 choose 2}+{n+1 choose 3}6+{n+1 choose 4}6$$

I like this visual proof, due to Man-Keung Siu. It appeared in the March 1984 issue of Mathematics Magazine.

See also two more proofs (as well as this one) in Roger Nelson's Proofs Without Words: Exercises in Visual Thinking.

Yet another proof(!)

Notice that $(k+1)^3 - k^3 = 3k^2 + 3k + 1$ and hence

$$(n+1)^3 = sum_{k=0}^n left[ (k+1)^3 - k^3right] = 3sum_{k=0}^n k^2 + 3sum_{k=0}^n k + sum_{k=0}^n 1$$

which gives you

$$begin{align}
sum_{k=1}^n k^2
& = frac{1}{3}(n+1)^3 - frac{1}{2}n(n+1) - frac{1}{3}(n+1) \
& = frac{1}{6}(n+1) left[ 2(n+1)^2 - 3n - 2right] \
& = frac{1}{6}(n+1)(2n^2 +n) \
& = frac{1}{6}n(n+1)(2n+1)
end{align}$$

Proof (by induction)

Basis: Check it for n = 1 (it works out).

Induction: Assume the result is true for a given value of $n$. That is, assume
$$
sum_{k = 1}^n k^2 = frac{n(n+1)(2n+1)}{6}.
$$
Try to show that the result holds for $n+1$.
$$
begin{align*}
sum_{k = 1}^{n+1} k^2 &= (n+1)^2 + sum_{k=1}^n k^2\
&= (n+1)^2 + frac{n(n+1)(2n+1)}{6}\
&= frac{6(n+1)^2 + n(n+1)(2n+1)}{6}\
&= frac{(n+1)(n+1+1)(2(n+1)+1)}{6}.
end{align*}
$$

To verify the identity, note $rm:sum_{k=1}^n: k^2 = f(n) iff f(n+1) - f(n) = (n+1)^2:$ and $rm: f(1) = 1:. $ But it's rote polynomial arithmetic to check that the RHS polynomial satisfies this recurrence.

To discover the identity, notice that any polynomial solution of the above recurrence has degree at most $3$. Hence it's easy to find the polynomial solution by substituting a cubic polynomial with undetermined coefficients.

Generally one can give a formula for sums of power using Bernoulli polynomials (motivated by discrete analogs of integrals of powers). The general theory becomes much clearer when one studies finite difference calculus and umbral calculus.

I think it's useful to report here another proof that I have posted on Mathoverflow.

Write down numbers in an equilateral triangle as follows:

    1

   2 2    

  3 3 3

 4 4 4 4

Now, clearly the sum of the numbers in the triangle is $Q_n:=1^2+2^2+dots+n^2$. On the other hand, if you superimpose three such triangles rotated by $120^circ$ each, like these ones

    1          4          4 

   2 2        3 4        4 3

  3 3 3      2 3 4      4 3 2

 4 4 4 4    1 2 3 4    4 3 2 1

then the sum of the numbers in each position equals $2n+1$. Therefore, you can double-count $3Q_n=frac{n(n+1)}{2}(2n+1)$. $square$

The proof is not mine and I do not claim otherwise. I first heard it from János Pataki. It is similar (but simpler) to the proof appearing on Wikipedia as I am writing this.

How to prove formally that all positions sum to $2n+1$? Easy induction: moving down-left or down-right from the topmost number does not alter the sum, since one of the three summand increases and one decreases. This is a discrete analogue of the Euclidean geometry theorem "given a point $P$ in an equilateral triangle $ABC$, the sum of its three distances from the sides is constant" (proof: sum the areas of $APB,BPC,CPA$), which you can mention as well.

A probabilistic method that I learned from Jim Pitman's book Probability (exercise 3.3.10) is as follows. Let $X$ be uniformly distributed on the set ${ 1, 2, ldots, n }$. Then
$$ E(X^3) = (1^3 + 2^3 + ldots + n^3)/n $$
and
$$ E((X+1)^3) = (2^3 + 3^3 + ldots +(n+1)^3)/n. $$
Subtracting the first of these from the second we get
$$ E((X+1)^3 - X^3) = ((n+1)^3 - 1)/n $$
and we can simplify both sides a bit to get
$$ E(3X^2 + 3X + 1) = n^2 + 3n + 3.$$
By linearity of expectation we can expand the left-hand side to get
$$ 3 E(X^2) + 3 E(X) + 1 = n^2 + 3n + 3. $$

Now $E(X) = (1+2+ldots+n)/n = (n+1)/2$. Substituting this in and solving for $E(X^2)$ gives

$$ E(X^2) = {(n+1)(2n+1) over 6} $$
but of course $E(X^2) = (1^2+2^2+cdots +n^2)/n$.

Similarly, we can derive for each $k$
$$ sum_{j=0}^{k-1} {k choose j} E(X^j) = sum_{l=1}^k {k choose l} n^{l-1} $$
and so if we know $E(X^0), ldots, E(X^{k-2})$ we can solve for $E(X^{k-1})$. So this method generalizes to higher moments as well.

Sums of polynomials can be done completely mechanically (no insight required, just turn the handle!) using the discrete calculus. Bill Dubuque mentions this in his answer, but I think it's nice to see a worked example.

Represent $k^2$ in terms of falling powers (easy by inspection in this case, but you can use Stirling subset numbers to convert): $$ k^2 = k^{underline 2} + k^{underline 1}$$

Sums of falling powers are easy, just like integration of ordinary powers, except for the treatment of limits: $$ sum_{k=1}^n k^{underline 2} + k^{underline 1} = bigg({1over 3}k^{underline 3} + {1over 2}k^{underline 2}bigg) bigg|^{n+1}_0$$

And then convert back into ordinary powers (by expansion, or using signed Stirling cycle numbers): $$ {1over 3}((n+1)^3 - 3(n+1)^2 + 2(n+1)) + {1over 2}((n+1)^2 - (n+1))$$

And then you can rearrange to get the answer you want.

This is a method that I learned from Polya's Mathematics and Plausible Reasoning: Let $s(n) = 1 + 2 + cdots + n$ and let $t(n) = 1^2 + 2^2 + cdots + n^2$. Make a small table as follows:

   n = 1 2  3  4  5

t(n) = 1 5 14 30 55

s(n) = 1 3  6 10 15

Note the ratio $r(n) = t(n)/s(n)$ for sucessive values of $n$:

R(1) = 1 = 3/3

R(2) = 5/3

R(3) = 14/6 = 7/3

R(4) = 30/10 = 3 = 9/3

R(5) = 55/15 = 11/3

Based on the pattern it seems that $r(n) = (2n+1)/3$ (and in fact it is: just prove it by induction). It follows that $t(n) = r(n)s(n)$. Now use the fact that $s(n) = n(n+1)/2$.

A natural approach for this kind of problems when you don't know the result is to proceed as follows :

We may want to write the sum $sum_{k=1}^n k^2$ as a telescopic sum, so we will try to find a polynomial of degree 3 ( why ? ) $P$ so that $Pleft( k+1 right) - Pleft(kright)=k^2$. Let $Pleft( x right) = ax^3+bx^2+cx$ for all reals $x$, then our constraint becomes :

$k^2= aleft( left(k+1right)^3 - k^3 right) + bleft( left(k+1right)^2 - k^2 right) + c$

Which after expanding and rearranging becomes :

$k^2 = 3ak^2 + left( 3a+2b right)k + a+b+c$

But we know that two polynomials are equal iff their coefficients are equal too, so we just need to solve this system :

$left{ begin{aligned}
a &= frac{1}{3} \
3a+2b &= 0 \
a+b+c &= 0
end{aligned} right.$

Which gives us $left( a,b,c right) = left( frac{1}{3}, frac{-1}{2}, frac{1}{6} right)$

And Voilà, we just found the coefficients of our polynomial ! Now we just have to evaluate our telescopic sum :

$sum_{k=1}^n k^2 = sum_{k=1}^n Pleft( k+1 right) - Pleft(kright) = Pleft(n+1right)-underbrace{Pleft(1right)}_{=0}$

$sum_{k=1}^n k^2 = frac{1}{3}left(n+1right)^3 - frac{1}{2}left(n+1right)^2+frac{1}{6}left(n+1right)$

$sum_{k=1}^n k^2 = frac{1}{6}left(n+1right)left( 2 left(n+1right)^2 - 3 left(n+1 right) + 1 right)$

$sum_{k=1}^n k^2 = frac{1}{6}left(n+1right)left( 2n^2+n right)$

$sum_{k=1}^n k^2 = frac{1}{6}nleft(n+1right)left(2n+1right)$

Which completes the proof :-)

The standard method is induction and you should look it up as it is a popular second example (first is $sum n$)

Another argument is use: $$24n^2 +2= (2n+1)^3-(2n-1)^3$$ and get a telescoping sum.

i.e $$24sum_1^n k^2 +2n = sum_1^n (2k+1)^3-sum_1^n (2k-1)^3$$ $$24sum_1^n k^2 +2n = (2n+1)^3-1$$
$$24sum_1^n k^2 =8 n^3+12 n^2+4 n$$ $$24sum_1^n k^2 =4 n (n+1) (2 n+1)$$ $$sum_1^n k^2 = frac{n (n+1) (2 n+1)}{6}$$

Proof 1. (Exercise 2.5.1 in Dias Agudo, Cândido da Silva, Matemáticas Gerais III). Let $S:=sum_{k=1}^{n}k^{2}$. Consider $(1+a)^{3}=1+3a+3a^{2}+a^{3}$ and sum
$(1+a)^{3}$ for $a=1,2,ldots ,n$:

$$begin{eqnarray*}
(1+1)^{3} &=&1+3cdot 1+3cdot 1^{2}+1^{3} \
(1+2)^{3} &=&1+3cdot 2+3cdot 2^{2}+2^{3} \
(1+3)^{3} &=&1+3cdot 3+3cdot 3^{2}+3^{3} \
&&cdots \
(1+n)^{3} &=&1+3cdot n+3cdot n^{2}+n^{3}
end{eqnarray*}$$

The term $(1+1)^3$ on the LHs of the 1st sum cancels the term $2^3$ on the RHS of the 2nd, $(1+2)^3$, the $3^3$, $(1+3)^4$, the $4^3$, ..., and $(1+n-1)^3$ cancels $n^3$. Hence

$$(1+n)^{3}=n+3left( 1+2+ldots +nright) +3S+1$$

and

$$S=frac{n(n+1)(2n+1)}{6},$$

because $1+2+ldots +n=dfrac{nleft( n+1right) }{2}$.

Proof 2. (Exercise 1.42 in Balakrishnan, Combinatorics, Schaum's Outline of Combinatorics). From

$$binom{k}{1}+2binom{k}{2}=k+2frac{kleft( k-1right) }{2}=k^{2},$$

we get

$$begin{eqnarray*}
S &:&=sum_{k=1}^{n}k^{2}=sum_{k=1}^{n}binom{k}{1}+2binom{k}{2}
=sum_{k=1}^{n}binom{k}{1}+2sum_{k=1}^{n}binom{k}{2} \
&=&binom{n+1}{2}+2binom{n+1}{3} \
&=&frac{nleft( n+1right) left( 2n+1right) }{6}.
end{eqnarray*}$$

Here's one using the Pertubation Method I learnt in Concrete Mathematics:
$$S_n = sum_{0leq jleq n}j^3$$.
$$S_n+(n+1)^3=sum_{0leq jleq n+1}j^3$$
$$S_n+(n+1)^3=0+sum_{1leq jleq n+1}j^3$$
Replacing $j$ by $j+1$ gives us
$$S_n+(n+1)^3=sum_{1leq j+1 leq n+1}(j+1)^3$$
Rewriting $1leq j+1leq n+1$ as $0leq jleq n$ and expanding$(j+1)^3$
$$S_n+(n+1)^3=sum_{0leq j leq n}(j^3+1+3j^2+3j)$$
By Associative law
$$S_n+(n+1)^3=sum_{0leq j leq n}j^3 +sum_{0leq jleq n}1 + 3sum_{0leq j leq n}j^2+3sum_{0leq jleq n}j$$
$S_n =sum_{0leq jleq n}j^3$, so it gets canceled.
Rewriting $sum_{0leq jleq n}1$ and $sum_{0leq jleq n}j$ as $(n+1)$ and $frac{n(n+1)}{2}$ respectively
$$(n+1)^3=(n+1)+frac{3n(n+1)}{2}+3sum_{0leq jleq n}j^2$$
$$3sum_{0leq jleq n}j^2=(n+1)^3-frac{3n(n+1)}{2} - (n+1)$$
$$3sum_{0leq jleq n}j^2=(n+1)(n^2+1+2n-frac{3n}{2}-1)$$
$$3sum_{0leq jleq n}j^2=frac{(n+1)(2n^2+n)}{2}$$
$$sum_{0leq jleq n}j^2=frac{n(n+1)(2n+1)}{6}$$
Using the same methods, one can get closed forms for even higher sums like $sum_{j=0}^{n}j^3$ by taking $S_n = sum_{0leq jleq n}j^4$ and using the binomial expansion for $(j+1)^4$

A combinatorial proof:

Let $S={1,2,dots,(n+1)},nge 2$ and $T={(x,y,z)|x.y,zin S,x< z,y< z}$.By counting the number of members of $T$ in $2$ different ways I will prove the formula.

$1$st way:

We will at first Choose $z$ form the set $S$.When $z$ is $1$ then there are no choices for $x,y$ so the no. of elements of $T$ with $z=0$ is zero.When $z=2$ the number of choices for $x$ is $1$ and so is for $y$(precisely $x=y=1$).When $z=3$ then $xin {1,2}$ and $yin {1,2}$ so total no. of choices equals $2^2$.In a similar manner when $z=k,(1le kle (n+1))$,no. of choices for $x$ equals $(k-1)$ and no. of choices for $y$ is also $(k-1)$. So total no . of elements of T with $z=k$ is $(k-1)^2$.

So we will get the total no. of elements of $T$ by summing $(k-1)^2$ up from $1 $ to $(n+1)$.Hence $$|T|=sum _{l=1}^{(n+1)}(l-1)^2=sum_{k=1}^{n}k^2$$

$2$nd way:

Among the elements of $T$ consisting of three numbers from the set $S$, there are elements in $x=y$ and elements in which $xne y$.

We can count the no. of elements in which $x=y$ by choosing two distinct nos. from $S$ and assigning $z$ with the lagest no. and $x,y$ with the smallest number. We can choose two distinct numbers from $S$ in $displaystyle binom{n+1}{2}$ ways, so the total no. elements having $x=y$ is $displaystyle binom{n+1}{2}$.

Now we have to count the number of elements in which $xne y$.This means that $x,y$ are dinstict and as they are less than $z$ this means that all the three are distinct. So we can count no. of such elements in $T$ in the following way.At first we will choose three elements from the set $S$ and assign the largest value to $z$ and assign the other two values to $x,y$. Now we can choose three numbers from the set $S$ in $displaystyle binom{n+1}{3}$.From each such three element we can get two elements of the set $T$(Assigning the largest to $z$ and then assigning any one of then to $x$ and the other to $y$). So no. of elements of $T$ having $xne y$ is $2displaystyle binom{n+1}{3}$

So by this method we have $|T|=displaystyle binom{n+1}{2}+2displaystyle binom{n+1}{3}$

On equating the result obtained from both the methods we have
$$sum_{k=1}^{n}k^2=displaystyle binom{n+1}{2}+2displaystyle binom{n+1}{3}=frac{n(n+1)(2n+1)}{6}$$

Note that this can easily be extended to find the sum of $p$th power of integers.($pin mathbb{N})$

Another simple proof could be as follows: note that each square
can be written as a sum of odd numbers:

$sum_{k=1}^n(2k-1)=n^2$.

(that can be easily shown)

When writing each square as a sum of odd numbers we get that

$S=sum_{k=1}^n k^2=1 + (1+3) + (1+3+5) + ... =sum_{k=1}^n(n-k+1)(2k-1)=$

$=(2n+3)sum_{k=1}^n k -(n+1)sum_{k=1}^n 1 -2S$.

Therefore

$3S=frac{(2n+3)n(n+1)}{2}-n(n+1)=frac{(2n+1)n(n+1)}{2}$.

$begin{aligned} & hspace{0.5in} begin{aligned}displaystyle sum_{1 le k le n}k^2 & = sum_{1 le k le n}~sum_{1 le r le k}r =sum_{1 le r le n}~sum_{r le k le n}r \& = sum_{1 le r le n}~sum_{1 le k le n}r-sum_{1 le r le n}~sum_{1 le k le r-1}r \& = nsum_{1 le r le n}r-frac{1}{2}sum_{1 le r le n}r(r-1) \& =frac{1}{2}n^2(n+1)-frac{1}{2}sum_{1 le r le n}r^2+frac{1}{2}sum_{1 le r le n}r \& =frac{1}{2}n^2(n+1)-frac{1}{2}sum_{1 le k le n}k^2+frac{1}{4}n(n+1) end{aligned} \& begin{aligned}impliesfrac{3}{2}sum_{1 le k le n}k^2 & = frac{1}{2}n^2(n+1)+frac{1}{4}n(n+1) \& = frac{1}{4}n(n+1)(2n+1) end{aligned}\& implies hspace{0.15in} displaystyle sum_{1 le k le n}k^2 = frac{1}{6}n(n+1)(2n+1).end{aligned}$

Yet another take. Start with the definition of falling factorial powers:

$begin{align}
x^{underline{r}}
= x (x - 1) dotsm (x - r + 1)
end{align}$

so that:

$begin{align}
Delta n^{underline{r}}
&= (n + 1)^{underline{r}} - n^{underline{r}} \
&= r n^{underline{r - 1}} \
sum_{0 le k < n} k^{underline{r}}
&= frac{n^{underline{r + 1}}}{r + 1}
end{align}$

(the last one is also easy to prove by induction, or otherwise).

Now note that:

$begin{align}
n^2
= n^{underline{2}} + n^{underline{1}}
end{align}$

so we can write:

$begin{align}
sum_{0 le k le n} k^2
&= sum_{0 le k < n + 1}
left( k^{underline{2}} + k^{underline{1}} right) \
&= frac{(n + 1)^{underline{3}}}{3}
+ frac{(n + 1)^{underline{2}}}{2} \
&= frac{(n + 1) n (n - 1)}{3} + frac{(n + 1) n}{2} \
&= frac{(2 n + 1) (n + 1) n}{6}
end{align}$

My contribution:
$$begin{align}
sum_{k=1}^n k^2&=frac 14sum_{k=1}^n(2k)^2\
&=frac 14sum_{k=1}^nbinom {2k}2+binom {2k+1}2\
&=frac 14sum_{k=2}^{2n+1} binom k2\
&=frac 14binom {2n+2}3\
&=frac 14cdot frac {(2n+2)(2n+1)(2n)}{1cdot 2cdot 3}
=color{lightgrey}{frac{(n+1)(2n+1)n}{1cdot 2cdot 3}}\
&=frac 16n(n+1)(2n+1)quadblacksquare
end{align}$$

Another take, similar to Euler's (?) proof given by @leonbloy. We know that if:

$begin{align}
A(z)
&= sum_{n ge 0} a_n z^n
end{align}$

then (writing $mathtt{D}$ for the derivative):

$begin{align}
z mathtt{D} A(z)
&= sum_{n ge 0} n a_n z^n
end{align}$

Also:

$begin{align}
frac{A(z)}{1 - z}
&= sum_{n ge 0} left( sum_{0 le k le n} a_n right) z^n \
frac{1}{1 - z}
&= sum_{n ge 0} z^n
end{align}$

This you can repeat and combine. In our case, we get that:

$begin{align}
sum_{n ge 0} n^2
&= (z mathtt{D})^2 frac{1}{1 - z} \
&= frac{z + z^2}{(1 - z)^3} \
sum_{n ge 0} left( sum_{0 le k le n} k^2 right) z^n
&= frac{z + z^2}{(1 - z)^4}
end{align}$

We are interested in the coefficient of $z^n$:

$begin{align}
[z^n] frac{z + z^2}{(1 - z)^4}
&= [z^n] frac{z}{(1 - z)^4} + [z^n] frac{z^2}{(1 - z)^4} \
&= [z^{n - 1}] (1 - z)^{-4} + [z^{n - 2}] (1 - z)^{-4} \
&= (-1)^{n - 1} binom{-4}{n - 1} + (-1)^{n - 2} binom{-4}{n - 2} \
&= binom{n - 1 + 4 - 1}{4 - 1} + binom{n - 2 + 4 - 1}{4 - 1} \
&= binom{n + 2}{3} + binom{n + 1}{3} \
&= frac{(n + 2) (n + 1) n}{3!} + frac{(n + 1) n (n - 1)}{3!} \
&= frac{(2 n + 1) (n + 1) n}{6}
end{align}$

This approach is less messy than the cited one (no horrible derivatives and then l'Hôpital thrice).

Another method:

$$begin{align}
sum_{k=1}^nk^2&=frac 12sum_{k=1}^n k(k+1)+(k-1)k\
&=frac 12 left[sum_{k=1}^n k(k+1)+sum_{k=1}^{n-1}k(k+1)right]\
&=color{lightgrey}{frac 12} left[ color{lightgrey}2sum_{k=1}^n binom {k+1}2+color{lightgrey}2sum_{k=1}^{n-1}binom {k+1}2right]\
&=binom {n+2}3+binom {n+1}3\
&=frac{color{blue}{(n+2)}(n+1)n}6+frac{(n+1)ncolor{blue}{(n-1)}}6\
&=frac {n(n+1)color{blue}{(2n+1)}}6quadblacksquare
end{align}$$

A High School Proof:
$$S_n = 1^2 + 2^2 + 3^2 +dots+ n^2$$

We know,

$$r^3 - 3r^2 + 3r - 1 = (r-1)^3 $$

$$r^3 - (r-1)^3 = 3r^2 - 3r + 1$$

When $r=1, 2, 3,dots, n$

$$1^3 - 0^3 = 3*1^2 - 3*1 + 1qquad(1)$$

$$2^3 - 1^3 = 3*2^2 - 3*2 + 1qquad(2)$$

$$3^3 - 2^3 = 3*3^2 - 3*3 + 1qquad(3)$$

$$dotsdotsdotsdotsdotsdotsdotsdotsdots$$

$$n^3 - (n-1)^3 = 3n^2 - 3n + 1qquad(n)$$

By summing all the equations (from 1 to n) we get,

$$n^3 - 0^3 = 3(1^2 + 2^2 + 3^2 +dots+ n^2) - 3(1+2+3+dots+n) + (1+1+1+dots)$$

$$n^3 = 3S_n - 3frac{n(n+1)}{2} + n$$

begin{align}
3S_n & = n^3 + frac{3n(n+1)}{2} - n\
& = frac{2n^3 + 3n^2 + 3n - 2n}{2}\
& = frac{2n^3 + 3n^2 + n}{2}\
& = frac{n(2n^2 + 3n + 1)}{2}\
& = frac{n(2n^2 + 2n + n + 1)}{2}\
& = frac{n{2n(n+1)+1(n+1)}}{2}\
& = frac{n(n+1)(2n+1)}{2}\
end{align}

$$therefore S_n = frac{n(n+1)(2n+1)}{6}$$

I will refer to Qiaochu's excellent answer here as proof that if we define

$$f(N):=sumlimits_{n=0}^N n^2$$

then $f$ is a polynomial of degree $3$.

It is easy to calculate the first few values of this sum. Namely,

$begin{align}
f(0) &= 0 \
f(1) &= 1 \
f(2) &= 5 \
f(3) &= 14
end{align}$

I claim that these four points are sufficient to uniquely determine $f$.

To wit, we have in general that

$$f(x)=sumlimits_{k=0}^3 c_kx^k$$

which when be combined with the four computed values above results in the following system of equations:

$$begin{pmatrix}
1 & 0 & 0 & 0 \
1 & 1 & 1 & 1 \
1 & 2 & 4 & 8 \
1 & 3 & 9 & 27
end{pmatrix}
begin{pmatrix}
c_0 \ c_1 \ c_2 \ c_3
end{pmatrix}
=
begin{pmatrix}
0 \ 1 \ 5 \ 14
end{pmatrix}$$

This matrix is a Vandermonde matrix which has a well-known determinant

$$begin{align}
det(V) &= (1-0)(2-0)(3-0)(2-1)(3-1)(3-2) \
&neq 0
end{align}$$

Because its determinant is nonzero, the matrix is invertible, and so we have

$$begin{pmatrix}
c_0 \ c_1 \ c_2 \ c_3
end{pmatrix}
=
V^{-1}cdotbegin{pmatrix}
0 \ 1 \ 5 \ 14
end{pmatrix}$$

from which $f(x)$ can be determined directly.

However, if you're like most people, inverting a $4times4$ matrix doesn't exactly tickle your fancy!

Luckily, now that we see that the interpolating cubic is unique, we could find it through the described matrix multiplication, but we would get to the same result if we proceeded a different route as well. This is where Lagrange polynomials come to the rescue.

Using the general formula, we have immediately that

$$begin{align}
f(x) &= 0cdot(dots)+1cdotfrac{x(x-2)(x-3)}{1(1-2)(1-3)}+5cdotfrac{x(x-1)(x-3)}{2(2-1)(2-3)}+14cdotfrac{x(x-1)(x-2)}{3(3-1)(3-2)} \
&= frac{1}{2}left(x^3-5x^2+6xright)-frac{5}{2}left(x^3-4x^2+3xright)+frac{14}{6}left(x^3-3x^2+2xright) \
&= frac{1}{6}left(2x^3+3x^2+xright) \
&= frac{1}{6}xleft(x+1right)left(2x+1right)
end{align}$$

You can generalize this approach to find expressions for $sum n^pquadforall pinmathbb{N}$.

Or, you know, there's always Faulhaber's formula.

Another way to prove this by induction goes as follows:

Base case: If $n=0$, then we have $0$ on the left hand side, and $0(0+1)(2(0)+1)/6=0$ on the right.

Induction step:

Consider the differences $L(j+1)-L(j)$, and $R(j+1)-R(j)$ where $L(j)$ indicates that we have $j$ for $n$ on the left hand side. Well, $L(j+1)-L(j)=(j+1)^2$, and
$$R(j+1)-R(j)=frac{(j+1)((j+1)+1))(2(j+1)+1)}{6} - frac{j(j+1)(2j+1)}{6}$$ which simplifies to $(j+1)^2$ also. So, the rates of change on both sides equal each other, and thus the induction step follows.

Here, we present an approach that uses only the sum of the arithmetic progression $sum_{k=1}^nk=frac12n(n+1)$.

Here, we note that $sum_{j=1}^k(1)=k$. Then, we can write

$$begin{align}
sum_{k=1}^nk^2&=sum_{k=1}^nksum_{j=1}^k(1)\\
&=sum_{j=1}^nsum_{k=j}^n,k\\
&=frac12sum_{j=1}^n(n+1-j)(j+n)\\
&=frac12sum_{j=1}^nleft(n(n+1)+j-j^2right)\\
&=frac12n^2(n+1)+frac14n(n+1)-frac12sum_{j=1}^nj^2\\
frac32sum_{k=1}^nk^2&=frac{(2n+1)n(n+1)}{4}\\
sum_{k=1}^nk^2&=frac{(2n+1)n(n+1)}{6}
end{align}$$

For proof by induction; these are the $color{green}{mathrm{three}}$ steps to carry out:

Step 1: Basis Case: For $i=1$: $$sum^{i=k}_{i=1} i^2=frac{1(1+1)(2times 1+1)}{6}= frac{2times 3}{6}=1$$ So statement holds for $i=1$.

Step 2: Inductive Assumption: Assume statement is true for $i=k$:

$$sum^{i=k}_{i=1} i^2=frac{k(k+1)(2k+1)}{6} $$

Step 3: Prove Statement holds for $i=k+1$. You need to prove that for $i=k+1$: $$sum^{i=k+1}_{i=1} i^2=color{blue}{frac{(k+1)(k+2)(2k+3)}{6}}$$

To do this you cannot use: $$sum^{i=k}_{i=n} i^2=color{red}{frac{n(n+1)(2n+1)}{6}} $$ as this is what you are trying to prove.

So what you do instead is notice that:
$$sum^{i=k+1}_{i=1} i^2= underbrace{frac{k(k+1)(2k+1)}{6}}_{text{sum of k terms}} + underbrace{(k+1)^2}_{text{(k+1)th term}}$$
$$sum^{i=k+1}_{i=1} i^2= frac{k(k+1)(2k+1)}{6}+frac{6(k+1)^2}{6}$$
$$sum^{i=k+1}_{i=1} i^2= frac{(k+1)left(k(2k+1)+6(k+1)right)}{6}$$
$$sum^{i=k+1}_{i=1} i^2= frac{(k+1)(2k^2+color{green}{7k}+6)}{6}=frac{(k+1)(2k^2+color{green}{4k+3k}+6)}{6}=frac{(k+1)left(2k(k+2)+3(k+2)right)}{6}=color{blue}{frac{(k+1)(k+2)(2k+3)}{6}}quad forall space k in mathbb{N}$$

Which is the relation we set out to prove. So the method is to substitute $i=k+1$ into the formula you are trying to prove and then use the inductive assumption to recover the $color{blue}{mathrm{blue}}$ equation at the end.

Here is a sketch with calculus.

$f(x) = x^2$ is a strictly monotonously increasing function. The sum is an upper Riemann sum of width 1 of this function. But if we shift it down one step we get a lower Riemann sum of width 1. The difference in integral is $$int_{n-1}^n f(x)dx - int_0^1 f(x)dx = n^3/3 - (n-3)^3/3 - 1 = n^2-3n-4$$

This is not quite $n^2$, but we can now turn the question into asking "if we had an offset in the integration, say $+alpha$ on each "step", what $alpha$ should we select to get as close to $n^2$ as possible?"

This is not a full answer, but intends on helping how to think to try and solve problems.

Firstly, calculate the sum $S=sum_{k=1}^n k(k+1)$ which is:

$S=1cdot 2+2cdot 3+ cdots +n(n+1)$,

multiplying $S$ by 3
we get:

$3S=1cdot 2cdot 3+2cdot 3cdot (4-1)+3cdot 4cdot (5-2)+ cdots +ncdot (n+1)cdot (n+2-(n-1))$

$3S=1 · 2 · 3 + 2 · 3 · 4 − 1 · 2 · 3 + 3 · 4 · 5 − 2 · 3 · 4 + · · ·
+ n(n + 1)(n + 2) − (n − 1)n(n + 1)$

This telescoping series collapses to yield:

$$3S=n(n+1)(n+2)$$

$$S=frac{n(n+1)(n+2)}{3}$$

On the other side we have:

begin{alignat*}{2}
&sum_{k=1}^n k(k+1)&&=sum_{k=1}^n k^2+k \
& &&=sum_{k=1}^n k^2+sum_{k=1}^n k \
&frac{n(n+1)(n+2)}{3}&&=sum_{k=1}^n k^2+frac{n(n+1)}{2} \
&sum_{k=1}^n k^2&&=frac{n(n+1)(n+2)}{3}-frac{n(n+1)}{2} \
&sum_{k=1}^n k^2&&=frac{n(n+1)(2n+1)}{6}
end{alignat*}

Another method by using telescoping sum :-
We know $(a+b)^3-a^3-b^3=3ab(a+b)$ , take

$a=k-1 , b=2$ , then $a+b=k+1$ and $(k+1)^3-(k-1)^3-2^3=6(k-1)(k+1)=6k^2-6$ ,

hence $(k+1)^3-(k-1)^3-8+6=(k+1)^3-k^3+k^3-(k-1)^3-2=6k^2$ , taking sum over $k$

from $1$ to $n$ we get , $sum_{k=1}^n [(k+1)^3-k^3] + sum_{k=1}^n [k^3-(k-1)^3] -sum_{k=1}^n 2 = 6 sum_{k=1}^nk^2$ , the first

sum on the left hand side is telescoping resulting in $(n+1)^3-1$ , the second sum is also telescoping resulting in $n^3-(1-1)^3=n^3$ , and the third sum is simply $2n$ , hence

$6 sum_{k=1}^nk^2=(n+1)^3-1+n^3-2n=2n^3+3n^2+n=n(n+1)(2n+1)$

If you know that the nth square is the sum of the first n odd numbers, you can rewrite each square in the above sum in that way and do a little bit of rearranging to get the desired identity. In general, knowing the value of (n + 1)^k - n^k allows you to write (n+1)^k as a telescoping sum of a polynomial of degree k-1, running over the values 1 through n. If you know the values of the sums of consecutive powers up to k-1, this allows you to find the sum of consecutive kth powers by substituting the polynomial sum for each kth power.