Bespoke, 287 bytes

lower=A,upper=B
to find some way to apply squaring with sums to integers and find the solution
calculate with A+B,do squaring,times3to soon consider
multiply by some variants lessened from A-B
a division with dozenths completes solution
naturally for it,I do formulas expressing answer N

A port of the Aheui answer by @cobaltp.

Vyxal 3, 7 bytes

↯2℈¨Π∑d

-1 byte by emanresu A

Vyxal It Online!

↯2℈¨Π∑d­⁡​‎‎⁡⁠⁡‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁢‏⁠‎⁡⁠⁣‏‏​⁡⁠⁡‌⁣​‎‎⁡⁠⁤‏⁠‎⁡⁠⁢⁡‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁢⁢‏‏​⁡⁠⁡‌⁢⁡​‎‎⁡⁠⁢⁣‏‏​⁡⁠⁡‌­
↯        # ‎⁡inclusive range implicit input
 2℈      # ‎⁢all 2-wide combinations without replacement
   ¨Π    # ‎⁣product of each
     ∑   # ‎⁤sum
      d  # ‎⁢⁡double
💎

Created with the help of Luminespire.

<script type="vyxal3">
↯2℈¨Π∑d
</script>
<script>
    args=[["5","9"],["91","123"],["1","10"]]
</script>
<script src="https://themoonisacheese.github.io/snippeterpreter/snippet.js" type="module"/>

APL(NARS), 21 chars

(+/k*2)-⍨2*⍨+/k←↑../⎕

It would traslate the 2 numbers in the range, assign that range to k, make the sums, exp, and difference.

  (+/k*2)-⍨2*⍨+/k←↑../⎕
⎕:
  5 9
970
  (+/k*2)-⍨2*⍨+/k←↑../⎕
⎕:
  91,123
12087152
  (+/k*2)-⍨2*⍨+/k←↑../⎕
⎕:
  1 10
2640

A more readable function would be this:

{k←↑../⍵⋄(+/k*2)-⍨(+/k)*2}

but is a little more long (26 chars).

Desmos, 40 36 bytes _{(thx to emanresu A's suggestion)}

A=[x...y]
f(x,y)=A.total^2-A^2.total

Try this online!

Swift 6, 51 44 bytes

{{($1-1)*$1*(3*$0*$0+$1-2)/12}($0+$1,$0-$1)}

Try it on SwiftFiddle!

Haskell (38 34 bytes)

Shaved 4 bytes thanks to @Whea Wizard

l#h=sum[l..h]^2-sum[x*x|x<-[l..h]]

Perl 5 -MList::Util=sum -p, 35 bytes

$_=(sum@;=$_..<>)**2-sum map$_**2,@

Try it online!

J-uby, 32 bytes

:!~|:-%[:sum|g= ~:*.D,:*&g|:sum]

Attempt This Online!

(The space between = and ~ is necessary to keep the Ruby lexer from interpreting them as the =~ operator.)

Explanation

:!~ | :- % [:sum | g = ~:*.D, :* & g | :sum]

:!~ |                                         # Range, then
      :- % [                ,              ]  # difference of...
            :sum | g = ~:*.D                  #   sum, then square (save the square function in g), and
                              :* & g | :sum   #   map with square, then sum

~:*.D is the equivalent of ->x{ x*x } and a byte shorter than ~:**&2, the equivalent of ->x{ x**2 }. The ~ in the former case duplicates the given argument as many times as the function takes, and D forces it to behave dyadically. (The ~ in the latter case reverses the order of the arguments, so binding 2 to ~:** produces x**2 rather than 2**x—the ~ operator is very overloaded in J-uby.)

Tcl, 75 76 80 bytes

proc S l\ u {while \$l<=$u {incr Q $l
append R -$l*$l
incr l}
expr $Q**2$R}

Try it online!

Arturo, 42 35 bytes

$=>[-*<=∑r∑map r:<=&..&'z->z*z]

Try it

Thunno 2, 8 bytes

IDS²s²S-

Attempt This Online!

Explanation

IDS²s²S-  # Implicit input
ID        # Inclusive range; duplicate
  S²      # Sum the list; square the sum
    s     # Swap so the range is back on top
     ²S   # Square each item; sum the list
       -  # Take the difference
          # Implicit output

Ly, 22 bytes

nnR&s&+:*>l[:*sp<l->]<

Try it online!

This is a straight forward interpretation of the rules...

nnR                     - get the start/end, generate the sequence of numbers
   &s                   - save the sequence...
     &+                 - sum it
       :*               - square the sum
         >l             - load the sequence on a new stack
           [        ]   - while the stack isn't empty
            :*          - square the number
              sp<l      - move the squared number to the other stack
                  ->    - subtract and go back to the second stack
                     <  - switch back to the accumulator stack
                        - the number on the stack prints automatically on exit

Pyt, 8 bytes

ŘĐĐƩ²↔·-

Try it online!

Takes two inputs on separate lines as max, min

Ř               implicit inputs; Řange
 ĐĐ             Đuplicate array on stack twice
   Ʃ²           Get square of sum
     ⇹·         Get sum of squares
       -        Subtract; implicit print

JavaScript, 49 bytes

Without recursion:

(b,e)=>eval('for(x=y=0;b<=e;y+=b*b++)x+=b;x*x-y')

Try it:

f=(b,e)=>eval('for(x=y=0;b<=e;y+=b*b++)x+=b;x*x-y')

console.log(f(5,9)); // 970
console.log(f(91,123)); // 12087152
console.log(f(1,10)); // 2640

Vyxal, 8 bytes

ṡ:²∑$∑²ε

Try it Online!

How?

ṡ:²∑$∑²ε
ṡ        # Inclusive range between (implicit) second input and (implicit) first input
 :       # Duplicate this range
  ²      # Square each of the duplicate
   ∑     # Summate
    $    # Swap so the other range is at the top
     ∑   # Summate
      ²  # Then square
       ε # Absolute difference between the two

Factor + math.unicode, 34 bytes

[ [a,b] dup Σ sq swap norm-sq - ]

Try it online!

• [a,b] Create a range from the input.
• dup Σ sq swap Square of sum.
• norm-sq Sum of squares.
• - Subtract.

K (oK), 24 bytes

{(t*t:+/s)-+/s*s:x_!1+y}

Try it online!

"Boring" but straightforward.

Python 2, 43 bytes

f=lambda a,b,s=0:b/a and 2*a*s+f(a+1,b,s+a)

Test it on Ideone.

How it works

Call the function defined in the specification \$g(a, b)\$. We have that

$$ \newcommand{\sumab}[2]{\sum_{a \le #1 \le b} #2 \:} \begin{align} g(a,b) & = \left( \sumab n n \right)^2 - \sumab n {n^2} \\ & = \sumab {i,j} {ij} - \sumab n {n^2} \\ & = \sumab {i<j} {ij} + \sumab {i=j} {ij} + \sumab {j<i} {ij} - \sumab n {n^2} \\ & = 2 \sumab {j<i} {ij} \\ & = 2 \sumab {i} {\sum_{a \le i < j} ij} \\ & = \sumab {i} {\left( 2i \sum_{a \le j < i} j \right)} \end{align} $$

Define the function \$f(x, y, s)\$ recursively as follows.

$$ f(x, y, s) = \begin{cases} 2xs + f(x+1, y, s+x) & \text{ if } x \le y \\ 0 & \text{ if} x > y \end{cases} $$

By applying the recurrence relation of \$f(a, b, 0)\$ a total of \$b - a\$ times, we can show that:

$$ \begin{align} f(a, b, 0) & = 2\cdot0 + f(a+1, b, a) \\ & = 2\cdot0 + 2\cdot(a+1)\cdot a + f(a+2, b, a+(a+1)) \\ & = 2\cdot0 + 2\cdot(a+1)\cdot a + 2\cdot(a+2)\cdot(a+(a+1)) + f(a+3, b, a+(a+1)+(a+2)) \\ & \vdots \\ & = \sumab {i} {\left( 2i \sum_{a \le j < i} j \right)} + f\left(b+1, b, \sumab j j\right) \\ & = \sumab {i} {\left( 2i \sum_{a \le j < i} j \right)} \\ & = g(a, b) \end{align} $$

This is the function f of the implementation. While b/a returns a non-zero integer, the code following and is executed, thus implementing the recursive definition of f.

Once b/a reaches 0, we have that b > a and the lambda returns False = 0, thus implementing the base case of the definition of f.

Aheui (esotope), 108 bytes(36 chars)

방빠방빠쌍다쌍상싸사타빠밤타빠받다파반다받상싸사빠따따다따따받밤따나망히

Try it online!

Takes two inputs from user(separated by line feed or whitespace) and print result of following function;

\$f(a,b)=((a-b)\cdot(a-b-1)\cdot(3\cdot(a+b)^2+a-b-2))/12\$

Pony, 73 bytes

fun f(m:U64,n:U64,s:U64=0):U64=>if m>n then 0 else(2*m*s)+f(m+1,n,m+s)end

sort of ungolfed

  fun apply(m: U64, n: U64, s: U64 = 0): U64 =>
    if m > n then 0
    else (2 * m * s) + apply(m + 1, n, m + s) end

Port of Neil's JavaScript solution.

Rust, 67 bytes

|b,e|{let mut t=0;let mut s=0;for i in b..=e{t+=i;s+=i*i;}(t*t)-s};

A closure that takes two integers and returns an integer.

Burlesque - 16 bytes

r@J++2?^j)S[++?-

r@                  range
  J                 dup
   ++               sum
     2?^            square it
        j           swap
         )S[        map (square)
            ++      sum
              ?-    subtract

Alternate versions:

r@JJ?*++j++S[?-ab
r@J++2?^jqS[ms?-
r@J++J?*jqS[ms?-
r@J++J?*jJ?*++?-
r@J)S[++j++S[j?-

^- and this is the art of golfing :(.

Tcl, 105 bytes

proc d a\ b {while {[incr a]<=$b} {set j $a
while \$j<=$b {incr p [expr 2*([incr j]-1)*($a-1)]}}
expr $p}

Try it online!

# Tcl, 112 bytes

proc d a\ b {set i [expr $a-1]
while \$i<$b {set j [incr i]
while \$j<$b {incr p [expr 2*[incr j]*$i]}}
expr $p}

Bash, 74 65 bytes

x=y=0;for i in `seq $1 $2`;{ x=$[x+i];y=$[y+i*i]; };echo $[x*x-y]

Try it online!

I'm no expert at bash golfing, but this works.
Managed to cut 9 bytes just by reading through the bash golfing tips thread

Sidef, 46 44 bytes

Translation of: my Factor solution

{|a,b|R=a..b;R.map{.sqr}.sum-R.sum.sqr->abs}

Try it online!

use it like x.run(4, 20) -> 38760

Previous

->(a,b,R=a..b){R.map{.sqr}.sum-R.sum.sqr->abs}

longer, explained

->                       # the next list is a function parameter list
(a, b, R = a .. b )      # R is a RangeNum over the a and b 
{                        # a Block
  return                 # return the whole expression
  R.map{ |n|             # map this Block over the range
    return n.sqr         # squared value
  }.sum                  # summed the range
  -                      # subtract
  R.sum.sqr              # sum of the range, squared
  ->abs                  # call Number.abs on everything on the left-hand side of ->
}                        # ok, we're done 

Wolfram Language (Mathematica), 18 bytes

Tr@#^2-#.#&@*Range

Try it online!

Composition of Tr@#^2-#.#& (square of sum minus sum of squares, where sum of squares is implemented via dot product) and Range (list of all integers between the two inputs). If you like seeing the operations in order, Range/*Tr@#^2-#.#& is equivalent.

As a bonus, also solves this challenge with no modifications, since Range can either take one argument (giving the range from 1 to the input) or two. (But it's not the most efficient solution possible there.)

Python 3, 66 65 bytes

lambda a,b:sum(x*y*2for x in range(a,b+1)for y in range(x+1,b+1))

Try it online!

with @Dennis formulae and w/o recursion, yet another Python solution (as I cannot comment @leaky-nun solution (not enough reputation)

PHP, 58 bytes

for([,$p,$q]=$argv;$p<=$q;$s+=$p++)$t+=$p*$p;echo$s*$s-$t;

Run with -nr or try it online.

C# (.NET Core), 62 bytes

(a,b)=>{int x=0,y=0;for(;a<=b;a++){x+=a;y+=a*a;}return x*x-y;}

Try it online!

Ungolfed:

(a, b) => {                 // takes in two integer inputs, delimited by a comma
    int x = 0, y = 0;       // initializes x and y
    for(; a <= b; a++)      // increment a until a equals b
    {
        x += a;             // add current a to x
        y += a * a;         // add square of current a to y
    }
    return (x * x) - y;     // return the difference of the square of the sums (x*x) and the sum of the squares (y)
}

Excel, 60 bytes

=(B1^2+B1-A1^2+A1)^2/4-(2*B1^3+3*B1^2+B1-2*A1^3+3*A1^2-A1)/6

Algebraic reshuffling of more verbose implementation: For inputs a and b:

(B1*(B1+1)-A1*(A1-1))/2          // (Sum of b) - (Sum of a-1)
(  )^2                           //  squared
B1*(B1+1)*(2*B1+1)/6             //  Sum of (b squared)
A1*(A1-1)*(2*A1-1)/6             //  Sum of (a-1 squared)

Red, 70 bytes

func[a b][s: t: 0 until[s: s + a t: a * a + t  b < a: a + 1]s * s - t]

Try it online!

C (gcc), 48 bytes

a,b;f(n,m){for(a=b=0;m/n;a+=n++)b+=n*n;n=a*a-b;}

Try it online!

Julia 0.6, 21 bytes

r->2triu(r*r',1)|>sum

Try it online!

Julia port of Luis Mendo's MATL answer. Comes out to the same bytecount as the other Julia answer by @gggg. Takes a Range similar to that answer, does a matrix multiply of it with its transpose to get all pairwise products in a matrix. Takes the upper triangular portion of that (the 1 argument is to denote the distance from the leading diagonal to start from), multiplies that by 2, and finally sums the values and returns that implicitly.

Python 2.7, 82 bytes 85 bytes

def a(l, u):
    c=0
    d=0
    for e in range(l, u+1):
        c+=e
        d+=e**2
    return c**2-d

Works in base 10.

Jelly, 6 bytes

rµÄḋḊḤ

Recent improvements to the Jelly language allow a compact implementation of the \$g\$ function from my Python answer.

Try it online!

How it works

rµÄḋḊḤ  Main link. Arguments: a, b (integers)


r       Range; yield R := [a, ..., b].
 µ      Begin a monadic chain with argument R.
  Ä     Accumulate; take the cumulative sum of R.
    Ḋ   Deque; yield [a+1, ..., b].
   ḋ    Take the dot product, ignoring the last term of the cumulative sum.
     Ḥ  Unhalve; double the result.

C (gcc), 67 bytes

f(n,m){double a,b;for(;n<=m;a+=n,b+=pow(n++,2));return pow(a,2)-b;}

Try it online!

Perl 6,  36 32  31 bytes

{([+] $_=@_[0]..@_[1])²-[+] $_»²}
{([+] $_=$^a..$^b)²-[+] $_»²}
{[+]($_=$^a..$^b)²-[+] $_»²}

Test it

Explanation:

{ # bare block with placeholder parameters $a and $b

  [+](# reduce with &infix:<+>
      # create a range, and store it in $_
      $_ = $^a .. $^b
  )²
  -
  [+] # reduce with &infix:<+>
    # square each element of $_ ( possibly in parallel )
    $_»²
}

Test:

#! /usr/bin/env perl6
use v6.c;
use Test;

my @tests = (
  (5,9) => 970,
  (91,123) => 12087152,
  (1,10) => 2640,
);

plan +@tests;

my &diff-sq-of-sum = {[+]($_=$^a..$^b)²-[+] $_»²}

for @tests -> $_ ( :key(@input), :value($expected) ) {
  is diff-sq-of-sum(|@input), $expected, .gist
}

1..3
ok 1 - (5 9) => 970
ok 2 - (91 123) => 12087152
ok 3 - (1 10) => 2640

Haskell, 36 bytes

m#n=sum[2*i*j|i<-[m..n],j<-[i+1..n]]

λ> m # n = sum [ 2*i*j | i <- [m..n], j <- [i+1..n] ]
λ> 5 # 9
970
λ> 91 # 123
12087152
λ> 1 # 10
2640

Note that

$$\left( \sum_{k=m}^n k \right)^2 - \sum_{k=m}^n k^2 = \cdots = \sum_{k_1=m}^n \sum_{k_2=m\\ k_2 \neq k_1}^n k_1 k_2 = \sum_{k_1=m}^n \sum_{k_2=k_1+1}^n 2 \,k_1 k_2$$

Japt, 10 bytes

õV
x²aUx ²

Try it

Julia 0.6, 21 bytes

r->sum(r)^2-sum(r.^2)

Returns an anonymous function that takes a Range, eg f(1:10), which then returns the result.

Try it online!

Whispers v2, 90 bytes

> Input
> Input
> 2
>> 1…2
>> L*3
>> Each 5 4
>> ∑6
>> ∑4
>> 8*3
>> 9-7
>> Output 10

Try it online!

R - 33 bytes

x=scan():scan();sum(x)^2-sum(x^2)

pass your lower limit to the first scan then your upper limit to the second scan

GeoGebra, 91 bytes

a(x)=(x²+x)/2
b(x)=x³/3+x²/2+x/6
c(x,y)=(a(y)-a(x))²
d(x,y)=b(y)-b(x)
c(x-1,y)-d(x-1,y)

Defines a function (probably e(x,y)) that computes the desired difference.
a(x) calculates the sum of natural numbers between 0 and x.
b(x) calculates the sum of the squares of the natural numbers between 0 and x.
c(x,y) first computes the sum of the natural numbers between x and y, then squares that sum.
d(x,y) calculates the sum of squares between b(x) and b(y).
The last line defines a multi-variable function that finishes the calculation. The function is automatically assigned a name, saving a few bytes.

CoffeeScript, 49 Bytes

f=(n,m,s=0)->if n>m then 0else 2*n*s+f(n+1,m,n+s)

Translation of Neil's answer

Output:

f 91, 123 == 12087152 => true

Fith, 52 bytes

{ 1 + range dup sum 2 pow swap { 2 pow } map sum - }

This is an anonymous function that takes the two numbers on the stack and leaves a single number.

Explanation:

{
    1 + range dup      2 ranges from a to b inclusive
    sum 2 pow          Sum one and square it
    swap               Bring a fresh range to the top
    { 2 pow } map sum  Square every element and sum the list
    -                  Subtract
}

Factor, 48 bytes

[ [a,b] [ [ sq ] map sum ] [ sum sq ] bi - abs ]

An anonymous function.

[ 
  [a,b] ! a range from a to b 
  [ 
    [ sq ] map sum ! anonymous function: map sq over the range and sum the result 
  ] 
  [ sum sq ] ! the same thing, in reverse order
  bi - abs   ! apply both anon funcs to the range, subtract them and abs the result
]

JavaScript (ES6), 50 37 bytes

f=(n,m,s=0)=>n>m?0:2*n*s+f(n+1,m,n+s)

Now a port of @Dennis♦'s Python solution.

Perl 5.10, 41 bytes

map{$s+=$_,$r+=$_**2}(<>..<>);say$s**2-$r

Input is given on 2 lines, for instance:

61
127

Try it here!

MATL, 9 bytes

&:&*XRssE

Try it online!

Explanation

&:   % Inclusive range between the two implicit inputs
&*   % Matrix of all pair-wise products
XR   % Upper triangular part of matrix, without the diagonal
ss   % Sum of all elements of the matrix
E    % Multiply by 2. Implicit display

Example

These are the partial results of each line for inputs 5 and 9:

1. &:

   5 6 7 8 9
   

2. &:&*

   25 30 35 40 45
   30 36 42 48 54
   35 42 49 56 63
   40 48 56 64 72
   45 54 63 72 81
   

3. &:&*XR

   0 30 35 40 45
   0  0 42 48 54
   0  0  0 56 63
   0  0  0  0 72
   0  0  0  0  0
   

4. &:&*XRss

   485
   

5. &:&*XRssE

   970
   

Python 2, 45 bytes

lambda a,b:(a+~b)*(a-b)*(3*(a+b)**2+a-b-2)/12

Closed form solution - not the shortest, but I thought it'd be worth posting anyway.

Explanation

Let p(n) be the nth square pyramidal number, and t(n) be the nth triangular number. Then, for n over the range a, ..., b:

• ∑n = t(b)-t(a-1), and
• ∑n² = p(b) - p(a-1)
• So (∑n)²-∑n² = (t(b)-t(a-1))² - (p(b) - p(a-1)).

This expression reduces to that in the code.

TI-BASIC, 19 bytes

Prompt N,M
randIntNoRep(N,M
sum(Ans)2-sum(Ans2

randIntNoRep gets the range (shuffled). The rest is pretty self explanatory.

JavaScript (ES6), 46 bytes

f=(x,y,s=0,p=0)=>x<=y?f(x+1,y,s+x,p+x*x):s*s-p

Actually, 10 bytes

u@x;;*@Σ²-

Try it online!

Explanation:

u@x;;*@Σ²-
u@x         range(a, b+1)
   ;;       two copies
     *      dot product with itself (equal to sum of squares)
      @Σ²   sum remaining copy, square sum
         -  subtract

Julia, 25 bytes

f(a,b,x=a:b)=sum(x)^2-x'x

This is a function that accepts two integers and returns a 1x1 integer array.

The approach is simple: Construct a UnitRange from the endpoints a and b and call it x, then sum x, square it, and subtract its norm, which is computed as transpose(x) * x.

Try it online! (includes all test cases)

Java, 84 77 characters, 84 77 bytes

7 bytes smaller due to Martin Ender and FryAmTheEggMan, thank you.

public int a(int b,int c){int e=0,f=0;for(;b<=c;e+=b,f+=b*b++);return e*e-f;}

Using the three test cases in the original post: http://ideone.com/q9MZSZ

Ungolfed:

public int g(int b, int c) {
    int e = 0, f = 0;
    for (; b <= c; e += b, f += b * b++);
    return e*e-f;
}

Process is fairly self-explanatory. I declared two variables to represent the square of the sums and the sum of the squares and repeatedly incremented them appropiately. Finally, I return the computed difference.

R, 23 bytes

sum(i:j)^2-sum((i:j)^2)

Does pretty much what you expect by taking i and j as inputs.

Pyke, 11 bytes

h1:Ds]MXXs-

Try it here!

h1:         - inclusive_range(input)
   Ds]      -     [^, sum(^)]
      MX    -    deep_map(^, <--**2)
         s  -   ^[1] = sum(^[1])
          - -  ^[0]-^[1]

Octave, 27 23 bytes

@(x,y)sum(z=x:y)^2-z*z'

Creates an anonymous function named ans which accepts two inputs: ans(lower, upper)

Online Demo

Explanation

Creates a row vector from x to y (inclusive) and stores it in z. We then sum all the elements using sum and square it (^2). To compute the sum of the squares, we perform matrix multplication between the row-vector and it's transpose. This will effectively square each element and sum up the result. We then subtract the two.

Matlab, 30 29 28 bytes

Using Suever's idea of norm gives us 2 bytes less

@(x,y)sum(x:y)^2-norm(x:y)^2

Old (simple) version:

@(x,y)sum(x:y)^2-sum((x:y).^2)

Actually, 11 bytes

u@x;♂²Σ@Σ²-

Try it online!

Python, 62 bytes

lambda a,b:sum(range(a,b+1))**2-sum(x*x for x in range(a,b+1))

Or, longer:

lambda a,b:sum(sum(x*y for y in range(a,b+1)if x-y)for x in range(a,b+1))

Ideone both!

Java, 95 bytes

public int d(int a,int b){int i=a,j,s=0;for(;i<=b;i++)for(j=a;j<=b;j++)s+=i==j?0:i*j;return s;}

Ideone it!

Haskell, 34 bytes

a#b=sum[a..b]^2-sum(map(^2)[a..b])

Usage example: 91 # 123 -> 12087152.

Nothing to explain.

Pyth, 11 bytes

s*M-F#^}FQ2

Try it online!

s*M-F#^}FQ2
       }FQ    Compute the range
      ^   2   Generate all pairs
   -F#        Remove those pairs who have identical elements
 *M           Product of all pairs
s             Sum.

Labyrinth, 28 24 bytes

?:?:}+=-:(:(#{:**+**#2/!

Try it online!

Explanation

Since loops tend to be expensive in Labyrinth, I figured the explicit formula should be shortest, as it can be expressed as linear code.

Cmd Explanation                 Stacks [ Main | Aux ]
?   Read M.                     [ M | ]
:   Duplicate.                  [ M M | ]
?   Read N.                     [ M M N | ]
:   Duplicate.                  [ M M N N | ]
}   Move copy to aux.           [ M M N | N ]
+   Add.                        [ M (M+N) | N ]
=   Swap tops of stacks.        [ M N | (M+N) ]
-   Subtract.                   [ (M-N) | (M+N) ]
:   Duplicate.                  [ (M-N) (M-N) | (M+N) ]
(   Decrement.                  [ (M-N) (M-N-1) | (M+N) ]
:   Duplicate.                  [ (M-N) (M-N-1) (M-N-1) | (M+N) ]
(   Decrement.                  [ (M-N) (M-N-1) (M-N-2) | (M+N) ]
#   Push stack depth.           [ (M-N) (M-N-1) (M-N-2) 3 | (M+N) ]
{   Pull (M+N) over from aux.   [ (M-N) (M-N-1) (M-N-2) 3 (M+N) | ]
:   Duplicate.                  [ (M-N) (M-N-1) (M-N-2) 3 (M+N) (M+N) | ]
*   Multiply.                   [ (M-N) (M-N-1) (M-N-2) 3 ((M+N)^2) | ]
*   Multiply.                   [ (M-N) (M-N-1) (M-N-2) (3*(M+N)^2) | ]
+   Add.                        [ (M-N) (M-N-1) (3*(M+N)^2 + M - N - 2) | ]
*   Multiply.                   [ (M-N) ((M-N-1)*(3*(M+N)^2 + M - N - 2)) | ]
*   Multiply.                   [ ((M-N)*(M-N-1)*(3*(M+N)^2 + M - N - 2)) | ]
#   Push stack depth.           [ ((M-N)*(M-N-1)*(3*(M+N)^2 + M - N - 2)) 1 | ]
2   Multiply by 10, add 2.      [ ((M-N)*(M-N-1)*(3*(M+N)^2 + M - N - 2)) 12 | ]
/   Divide.                     [ ((M-N)*(M-N-1)*(3*(M+N)^2 + M - N - 2)/12) | ]
!   Print.                      [ | ]

The instruction pointer then hits a dead end and has to turn around. When it now encounters / it attempts a division by zero (since the bottom of the stack is implicitly filled with zeros), which terminates the program.

J, 29 bytes

Port of Doorknob's Jelly answer.

[:(+/@(^&2)-~2^~+/)[}.[:i.1+]

Usage

>> f = [:(+/@(^&2)-~2^~+/)[}.[:i.1+]
>> 91 f 123x
<< 12087152

Where >> is STDIN, << is STDOUT, and x is for extended precision.

JavaScript (ES6), 67 bytes

a=>b=>([s=q=0,...Array(b-a)].map((_,i)=>q+=(s+=(n=i+a),n*n)),s*s-q)

Test Suite

f=a=>b=>([s=q=0,...Array(b-a)].map((_,i)=>q+=(s+=(n=i+a),n*n)),s*s-q)
e=s=>`${s} => ${eval(s[0])}` // template tag format for tests
console.log(e`f(5)(9)`)
console.log(e`f(91)(123)`)
console.log(e`f(1)(10)`)

Mathematica, 21 bytes

Tr[x=Range@##]^2-x.x&

An unnamed function taking two arguments and returning the difference. Usage:

Tr[x=Range@##]^2-x.x&[91, 123]
(* 12087152 *)

There's three small (and fairly standard) golfing tricks here:

• ## represents both arguments at once, so that we can use prefix notation for Range. Range@## is shorthand for Range[##] which expands to Range[a, b] and gives us an inclusive range as required.
• Tr is for trace but using it on a vector simply sums that vector, saving three bytes over Total.
• x.x is a dot product, saving four bytes over Tr[x^2].

05AB1E, 8 bytes

ŸDOnsnO-

Explained

ŸD       # range from a to b, duplicate
  On     # sum and square first range
    s    # swap top 2 elements
     nO  # square and sum 2nd range
       - # take difference

Try it online

Jelly, 9 8 bytes

rµS²_²S$

Try it online!

r         inclusive range from first input to second input
 µ        pass the range to a new monadic chain
  S       the sum
   ²      squared
    _     minus...
     ²S$  the squares summed

Thanks to FryAmTheEggman for a byte!

PowerShell v2+, 47 bytes

Two variations

param($n,$m)$n..$m|%{$o+=$_;$p+=$_*$_};$o*$o-$p

$args-join'..'|iex|%{$o+=$_;$p+=$_*$_};$o*$o-$p

In both cases we're generating a range with the .. operator, piping that to a loop |%{...}. Each iteration, we're accumulating $o and $p as either the sum or the sum-of-squares. We then calculate the square-of-sums with $o*$o and subtract $p. Output is left on the pipeline and printing is implicit.

MATL, 11 bytes

&:ts2^w2^s-

Try it online!

Explanation:

&:           #Create a range from the input
  t          #Duplicate it
   s2^       #Sum it and square it
      w      #swap the two ranges
       2^s   #Square it and sum it
          -  #Take the difference

APL, 23 20 bytes

-/+/¨2*⍨{(+/⍵)⍵}⎕..⎕

Works in NARS2000.

Brachylog, 24 bytes

:efL:{:2^.}a+S,L+:2^:S-.

Expects the 2 numbers in Input as a list, e.g. [91:123].

Explanation

:efL                     Find the list L of all integers in the range given in Input
    :{:2^.}a             Apply squaring to each element of that list
            +S,          Unify S with the sum of the elements of that list
               L+:2^     Sum the elements of L, then square the result
                    :S-. Unify the Output with that number minus S

CJam, 17 bytes

q~),>_:+2#\2f#:+-

Test it here.

Explanation

q~       e# Read and evaluate input, dumping M and N on the stack.
),       e# Increment, create range [0 1 ... N].
>        e# Discard first M elements, yielding [M M+1 ... N].
_        e# Duplicate.
:+2#     e# Sum and square.
\2f#:+   e# Swap with other copy. Square and sum.
-        e# Subtract.

Alternatively, one can just sum the products of all distinct pairs (basically multiplying out the square of the sum, and removing squares), but that's a byte longer:

q~),>2m*{)-},::*:+