Vyxal, 17 bytes

∆/J:\#*ε‛+=Y4IYf§

Try it Online!

I thought I'd answer this with modern Vyxal.

∆/                # Take the hypotenuse
  J               # And append it to the other two
    \#*           # For each, repeat # that many times
   :   ε          # Repeat ^ that many times
        ‛+=Y      # Interleave with "+="
            4IY   # Interleave with four spaces
               f§ # Flatten and transpose, padding with spaces

05AB1E, 38 34 bytes

nOtª©Å1'#×ζε®×R„= NĀèð.øý}»R„=+`.;

Takes the input as pair of integers (i.e. [3,4]).

Try it online or verify all test cases.

Explanation:

n            # Take the square of each value in the (implicit) input-list
 O           # Take the sum of that list
  t          # Take the square-root of that sum
   ª         # Append it to the (implicit) input-list
             #  i.e. [3,4] → [3,4,5.0]
    ©        # Store it in the register (without popping)
Å1           # Change each value to an inner list of that amount of 1s
  '#×       '# And then each 1 to a "#"
             #  → [["#","#","#"],["#","#","#","#"],["#","#","#","#","#"]]
     ζ       # Zip/transpose; swapping rows/columns, with space as default filler
             #  → [["#","#","#"],["#","#","#"],["#","#","#"],[" ","#","#"],[" "," ","#"]]
ε            # Map the inner lists to:
 ®           #  Push the list from the register
  ×          #  Repeat the characters that many times
             #   i.e. [" ","#","#"] and [3,4,5.0] → ["   ","####","#####"]
   R         #  Reverse this list
 „=          #  Push string "= "
    NĀ       #  Push the map-index truthified (0 remains 0; everything else becomes 1)
             #   i.e. 0 → 0
             #   i.e. 3 → 1
      è      #  Use it to index into the string
       ð.ø   #  Surround it with spaces
             #   → "=" → " = "
             #   → " " → "   "
          ý  #  Join the map-list together with this string as delimiter
             #   → "#####   ####      "
}»           # After the map, join everything by newlines
             #  → "##### = #### = ###\n#####   ####   ###\n#####   ####   ###\n#####   ####      \n#####             "
  R          # Reverse the string
             #  → "             #####\n      ####   #####\n###   ####   #####\n###   ####   #####\n### = #### = #####"
   „=+`.;    # And replace the first "=" with "+"
             #  → "             #####\n      ####   #####\n###   ####   #####\n###   ####   #####\n### + #### = #####"
             # (after which the the result is output implicitly)

C (gcc), 122 119 bytes

-3 thanks to ceilingcat

#define P(n);for(i=0;i*i++<n;puts(""))for(j=0;j*j++<n;)putchar(35);puts(
i,j;f(a,b){P(a*a)"+")P(b*b)"=")P(a*a+b*b)"");}

Try it online!

Husk, 36 bytes

↔T' `+₁√ṁ□¹`+C1" = "JC1" + "m₁
SRR'#

Try it online!

This can be made shorter if I can find a good way to add the symbols in.

Vyxal, L, 47 43 38 bytes (old codepage)

λ\#*wnř;:£f⍎÷ð\+ð?¥f⍎÷ð\=ð??Ķĵ∑√I¥f⍎÷W

This was intense.

Also, according to code-golf's statistics, only a small percentage of people who see my answers actually upvote them. So if you enjoy this answer, consider upvoting, it's free, and you can change your mind at any time (given you have ≥2k rep). Enjoy the answer.

Explained

λn\#*wnř;:£f⍎ðJ\+JðJ?¥f⍎JðJ\=JðJ??Ķĵ∑√I¥f⍎J

λ                                               # Start a lambda that:
 n\#*                                           #   Repeats the character "#" n times THEN
     w                                          #   Wraps that into a stack AND
      nř                                        #   Repeats the list n times
        ;                                       # Close lambda
         :£f                                    # Assign the lambda to variable 'f' leaving an extra copy
            ⍎                                   # Apply that lambda to the first input
             ðJ\+JðJ                            # Append " ", "+" and " " to the result from above
                    ?¥f⍎                        # Apply the lambda to the second input
                        J                       # Add it to the current list
                         ðJ\=Jð                 # Append " ", "=" and " " to the result from above
                               ??Ķ              # Pair the two inputs into a list
                                  ĵ∑            # Square each number and sum the list (a^2 + b^2)
                                    √I          # Square root that number and convert it to an integer
                                      ¥f⍎       # Apply the lambda to that result (c)
                                         J      # And add it to the big list
                                                # The L flag then auto vertical joins the list and outputs it

CJam, 49 bytes

" +   = "S/3/[q~_2$mh:H]_'#f*:a.*.\:+SH*f.e|zW%N*

Try it online in the CJam interpreter.

How it works

" +   = "S/3/ e# Split at spaces, then into chunks of length 3.
              e# This pushes [["" "+" ""] ["" "=" ""]].
[             e#
  q~          e# Read and interpret all input from STDIN.
  _2$         e# Copy both integers.
  mh          e# Calculate the hypotenuse of the triangle with those catheti.
  :H          e# Save the result in H.
]             e# Collect catheti and hypotenuse in an array.
_'#f*         e# Copy and replace each length with a string of that many hashes.
:a            e# Wrap each string in an array.
.*            e# Vectorized repetition. Turns strings into square arrays.
.\            e# Interleave with the string of operators.
:+            e# Concatenate to form an array of strings.
SH*           e# Push a string of spaces of length H.
f.e|          e# Mapped vectorized logical OR; pads all strings with spaces to
              e# length H.
zW%           e# Zip and reverse; rotates the array.
N*            e# Join the strings, separating by linefeeds.

C# (.NET Core), 221, 194 bytes

This feels way too long. This version just loops to construct the string.

EDIT: Ascii-Only with a nice -27 byte golf using the string constructor for serial char additions! Also, ty for pointing out I was using Math.Sqrt not System.Math.Sqrt. This has been adjusted!

(a,b)=>{int c=(int)System.Math.Sqrt(a*a+b*b),j=c;var s="";while(j>0)s+=new string(j>a?' ':'#',a)+(j>1?"   ":" + ")+new string(j>b?' ':'#',b)+(j-->1?"   ":" = ")+new string('#',c)+"\n";return s;}

Try it online!

PowerShell, 139 137 135 bytes

^{-2 thanks to ASCII-only
-2 thanks to Mazzy}

param($a,$b)($c=[math]::sqrt($a*$a+$b*$b))..1|%{(($m=" ","#")[$_-le$a]*$a)," +"[$_-eq1],($m[$_-le$b]*$b)," ="[$_-eq1],("#"*$c)-join" "}

Try it online!

Calculating $c hurt and there's probably a better way to conditionally swap between # and . Builds a list of chunks and joins them together while conditionally adding the signs.

Perl 6, 99 bytes

{$!=sqrt $^a²+$^b²;flip map({map {[' ','#'][$^d>$_]x$d,' =+ '.comb[!$_*++$ ]},$!,$b,$a},^$!)X"
"}

Try it online!

Anonymous code block that takes two numbers and returns the full string with a leading newline and three leading spaces and one trailing on each line.

If we can use other characters instead of #, then I can save a byte by replacing '#' with \*.

Charcoal, 24 bytes

⊞θ₂ΣＸθ²Ｆ =+«←←←ι←Ｇ↑←↓⊟θ#

Try it online! Link is to verbose version of code. Takes input as an array of two elements. Explanation:

⊞θ₂ΣＸθ²

Append the hypotenuse to the inputs.

Ｆ =+«

Loop over the characters that appear to the right of each square in reverse order.

←←←ι←

Print that character leftwards with spacing.

Ｇ↑←↓⊟θ#

Pop the last number from the array and print a square of #s of that size.

Japt, 28 bytes

Takes input as an array of integers.

pUx²¬)ËÆDç'#
í"+="¬ûR3)c ·z3

Try it

                    :Implicit input of array U=[a,b]
pUx²¬)ËÆDç'#
p                   :Push
 U ²                :  Square each element in U
  x                 :  Reduce by addition
    ¬               :  Square root
     )              :End push
      Ë             :Map each D
       Æ            :  Map the range [0,D)
        Dç'#        :    Repeat "#" D times
í"+="¬ûR3)c ·z3
í                   :Interleave
 "+="¬              :  Split the string "+=" to an array of characters
      û             :  Centre pad each
       R3           :    With newlines to length 3
         )          :End interleave
          c         :Flatten
            ·       :Join with newlines
             z3     :Rotate clockwise 270 degrees

APL (Dyalog Extended), 33 29 bytes^SBCS

-3 due to my extensions of Dyalog APL.

Anonymous prefix lambda:

{⊖⍕,' +=',⍪{⍵ ⍵⍴⍕#}¨⍵,√+/⍵*2}

Try it online!

{…} "dfn"; ⍵ is the argument (side lengths)

 ⍵*2 square

 +/ sum

 √ square-root

 ⍵, prepend argument

 {…}¨ apply the following anonymous lambda to each

  # root namespace

  ⍕ format as text

  ⍵ ⍵⍴ use argument twice to reshape into matrix with those dimensions.

 ⍪ make into column

 ' ++=', prepend these three characters to the three rows

 , ravel (combine rows into list)

 ⍕ format as text

 ⊖ flip upside-down

Lua5.2, 257 241 227 222 bytes

r=io.read
a=r"*n"b=r"*n"c=math.sqrt(a^2+b^2)d=a+b
w=io.write
for i=1,c do
for j=0,d+c+5 do
w((j>d+5 or(i>c-b and j>a+2 and j<d+3)or(i>c-a and j<a))and"#"or(i==c and(j==a+1 and"+"or(j==d+4 and"="or" "))or" "))end
w"\n"end

• Edit1: Simplified reading
• Edit2: Removed more whitespaces
• Edit3: aliases abstraction of io functions inspired by another answer

JavaScript ES6, 155 134 140 129 bytes

(n,m)=>eval("for(o='',q=(b,s)=>' #'[z<b|0].repeat(b)+(z?'   ':s),z=i=Math.hypot(n,m);z--;)o+=q(n,' + ')+q(m,' = ')+q(i,'')+`\n`")

I've rewritten this with for. Lots of golfing still...

If something isn't working, let me know. I'll fix it in the morning.

Tested on Safari Nightly

Ungolfed:

(n,m)=>
   Array(
     z=Math.hypot(n,m)
   ).fill()
   .map((l,i)=>
      (q=(j,s)=>
        (z-i<=j?'#':' ')
        .repeat(j)+
         (z-i-1?' ':s)
      )
      (n,`+`)+
      q(m,`=`)+
      q(z,'')
   ).join`
   `

Explanation:

(Not updated) but still accurate enough.

(n,m)=> // Function with two arguments n,m
   Array( // Create array of length...
    z=Math.hypot(n,m) // Get sqrt(n^2+m^2) and store in z
   ).fill() // Fill array so we can loop
   .map((l,i) => // Loop z times, take l, and i (index)
     (q=j=>( // Create function q with argument j
      z-i<=j? // If z-i is less than or equal to j...
        '#' // Use '#'
      : // OR
        ' ' // Use space
      ).repeat(j) // Repeat the character j times
     )(n) // Run with n
   + // Add to string
   ` ${ // Space
      (b=z-i-1)? // If this isn't the last line...
       ' ' // Return ' '
      : // Otherwise
       '+' // Plus
    } ${ // Space
      q(m) // run function q with arg m
    } ${ // Space
      b? // If b
       ' ' // Return space
      : // Otherwise
        '=' // '='
    }` + // Add to...
    '#'.repeat(z) // Repeat hashtag, z times
  ).join` // Join the new array with new lines
  `

DEMO

ES5 version Input must be valid sets of numbers:

function _taggedTemplateLiteral(e,t){return Object.freeze(Object.defineProperties(e,{raw:{value:Object.freeze(t)}}))}var _templateObject=_taggedTemplateLiteral(["\n"],["\n"]),t=function(e,t){return Array(z=Math.sqrt(e*e+t*t)).fill().map(function(r,n){return(q=function(e,t){return(z-n<=e?"#":" ").repeat(e)+(z-n-1?" ":t)})(e,"+")+q(t,"=")+q(z,"")}).join(_templateObject)};
// Demo
document.getElementById('go').onclick=function(){
  document.getElementById('output').innerHTML = t(+document.getElementById('input').value,
                                                 +document.getElementById('input2').value)
};

<div style="padding-left:5px;padding-right:5px;"><h2 style="font-family:sans-serif">Visually Explaining the Pythagorean Theorem</h2><div><div  style="background-color:#EFEFEF;border-radius:4px;padding:10px;"><input placeholder="Number 1" style="resize:none;border:1px solid #DDD;" id="input"><input placeholder="Number 2" style="resize:none;border:1px solid #DDD;" id="input2"><button id='go'>Run!</button></div><br><div style="background-color:#EFEFEF;border-radius:4px;padding:10px;"><span style="font-family:sans-serif;">Output:</span><br><pre id="output" style="background-color:#DEDEDE;padding:1em;border-radius:2px;overflow-x:auto;"></pre></div></div></div>

C, 176 bytes

C is not going to win this, but the fun is worth it.

#define A(x,y)for(j=x;j--;)putchar("# "[i+1>x]);printf(i?"   ":" "#y" ");
i;j;main(a,b,c){for(c=scanf("%d %d",&a,&b);a*a+b*b>c*c;c++);for(i=c;i--;puts("")){A(a,+)A(b,=)A(c,)}}

Pretty printed:

#define A(x,y)for(j=x;j--;)putchar("# "[i+1>x]);printf(i?"   ":" "#y" ");
i;j;
main(a,b,c)
{
    for(c=scanf("%d %d",&a,&b);a*a+b*b>c*c;c++);
    for(i=c;i--;puts(""))
    {
        A(a,+)
        A(b,=)
        A(c,)
    }
}

gcc enables us to pass third parameter to main (an array of environment variables), so we take advantage of it to use it for our purpose.

The

for(c=scanf("%d %d",&a,&b);a*a+b*b>c*c++;);

would be equivalent to

scanf("%d %d",&a,&b);
for(c=2;a*a+b*b>c*c++;);

because scanf returns the number of successfully scaned parameters.

Julia, 121 114 112 bytes

f(a,b)=for i=1:(c=isqrt(a^2+b^2)) g(x,t)=(i>c-x?"#":" ")^x*(i<c?"  ":t)" ";println(g(a," +")g(b," =")g(c,""))end

Ungolfed:

function f(a,b)
    # Compute the hypotenuse length
    c = isqrt(a^2 + b^2)

    # Write the lines in a loop
    for i = 1:c
        # Make a function for constructing the blocks
        g(x,t) = (i <= c - x ? " " : "#")^x * (i < c ? "  " : t) " "

        println(g(a," +") g(b," =") g(c,""))
    end
end

Fixed issue and saved 2 bytes thanks to Glen O.

Pyth, 35 32 31 30 bytes

j_.ts.i.imm*d\#d+Qs.aQ"+="mk4d

Try it online.

PHP, 178 170 168 bytes

Input is GET parameters x and y. Unfortunately I can't seem to golf those repeating strings.

<?php for(@$i=$z=hypot($x=$_GET[x],$y=$_GET[y]),@$s=str_repeat;$i;$i--)@print$s($i<=$x?~Ü:~ß,$x).(($l=$i==1)?~ßÔß:~ßßß).$s($i<=$y?~Ü:~ß,$y).($l?~ßÂß:~ßßß).$s(~Ü,$z).~õ;

• Saved 8 bytes by inverting all my strings and dropping the quotes.
• Saved 2 bytes by replacing the condition $i>0 with $i

Not sure why PHP doesn't like @echo so I had to sacrifice 1 byte with @print.

In case SE screws up the encoding, this is meant to be encoded in Windows-1252 (not UTF8).

Python 2, 134 100 bytes

a,b=input()
i=c=int(abs(a+b*1j))
while i:print"# "[i>a]*a," +"[i<2],"# "[i>b]*b," ="[i<2],"#"*c;i-=1

Try it online.

The program takes input as comma-separated integers, calculates the hypotenuse using Python's built-in complex numbers, then loops down from that value calculating and printing each line as it goes. The main golfing trick is using string indexing in place of conditionals to select #/+/= vs space.

Edit: The first version was a victim of some serious over-engineering--this one is both simpler and much shorter.

Pyth, 51 49 bytes

AQJs.aQLj*b]*b\#;j_MCm_.[d\ Jcj[yJb\=byHb\+byG))b

Expects input in the form [3,4].

Try it here

AQ - assigns input to G, H

Js.a,GH - calculates hypotenuse as J

Lj*b]*b\#; - defines y(b) as making a square of size b (elsewhere in the code, b means newline)

j_MCm_.[d\ Jcj[yJb\=byHb\+byG))b - Creates the squares, pads with spaces, and transposes

Saved two bytes thanks to Maltysen.

Ruby, 134

->a,b{c=((a**2+b**2)**0.5).round
c.times{|i|
d=i<c-1?'  ':'+='
puts (c-i>a ?' ':?#)*a+" #{d[0]}  #{(c-i>b ?' ':?#)*b} #{d[1]} "+?#*c}}

simple line by line approach.

Below in test program, with symbol changed to @ to help avoid confusting with the syntax #{....} ("string interpolation") used to insert expressions into a string. Each input should be given on a different line.

f=->a,b{c=((a**2+b**2)**0.5).round
c.times{|i|
d=i<c-1?'  ':'+='
puts (c-i>a ?' ':?@)*a+" #{d[0]}  #{(c-i>b ?' ':?@)*b} #{d[1]} "+?@*c}}

A=gets.to_i
B=gets.to_i
f.call(A,B)

CJam, 78 bytes

q~_2f#~+mQ+ee_2=~e>f{\~@1$-S*\'#*+_'#e=\a*_0=,S*@"+= "=1$,(S*\+1$a\a@a+++~}zN*

It first computes the hypotenuse (H), then, for each side (S), it builds an array of S lines made of: H-S spaces + S dashes. Finally, it transposes the matrix.

Demo