Wolfram Language (Mathematica), 104 93 63 bytes

g&[g=0;g+={{#2-g-{#2-g,c}~MinValue~x,c=-#<x<#}}&@@@#]/.x->0&

Try it online!

Input [{{radius, polynomial}...}], a list of pairs representing bowls with polynomials in \$x\$.

 is \[Piecewise]. An expression using If would have the same length but not evaluate correctly.

For decimal instead of symbolic output, use NMaxValue instead (or just call N on the result).

Builds a stack of bowls and evaluates the final one at \$x=0\$. The intermediate stacks make a nice visualization:

   g=0                                                        initial stack: height 0
                                                 @@@#         for each bowl:
       g+={{#2-g                    ,  -#<x<#}}                add newest bowl to stack
                 -{#2-g,c}~MinValue~x c=                        with offset
g&[                                                  ]/.x->0  height of stack at 0

Jelly, 54 53 bytes

J×$ÆrAƑƇ«⁹;⁹*€J{ḋ⁸ŻṀ
Œcz€0ḢṂç@I0;ⱮFƲƲ€ṚṁL’R€Ɗ;+Ṁ¥¥ƒ0Ṁ

Try it online!

A monadic link that takes as its argument the list of bowls from top to bottom in the format [[b1_radius, b1_coef1, ...], [b2_radius, b2_coef1, ...]] and returns the y position of the bottom of the top bowl.

Now correctly handles bowls which meet at places other than the minimal radius.

Explanation

Helper link: takes as left argument l the differences in coefficients of the polynomials representing the bowls from 1 upwards, and its right argument r the minimum radius; returns the maximum y value where the two bowls meet

  $                   | Following as a monad:
J                     | - Sequence from 1..<len(l)>
 ×                    | - Multiply (by l)
   Ær                 | Roots of polynomial
     AƑƇ              | Keep only those invariant when passed through absolute function (excludes negative, imaginary and complex numbers)
        «⁹            | Min of these filtered roots and r
          ;⁹          | Concatenate r to the list
            *€        | Each root/radius to the power of:
              J{      | - Sequence from 1..<len(l)>
                ḋ⁸    | Dot product with l
                  Ż   | Prepend zero
                   Ṁ  | Maximum

Main link, takes a bowl pile as its argument and returns y value of base of top bowl

Œc                               | Combinations length 2
  z€0                            | Transpose each using 0 as a filler
               Ʋ€                | For each one, do the following as a monad:
     Ḣ                           | - Take the head (the radii)     
      Ṃ                          | - Minimum
       ç@     Ʋ                  | - Call the helper link with this (min radius) as right argument and the following as left argument:
         I                       |   - Increments (difference between second and first polynomial for each coefficient)
          0;Ɱ                    |   - Prepend each with a zero (odd coefficients are all zero)
             F                   |   - Flatten
                 Ṛ               | Reverse
                  ṁ    Ɗ         | Mould to the following as a monad:
                   L             | Length
                    ’            | Decrease by 1
                     R€          | Range of each (e.g. [1],[1,2],[1,2,3],[1,2,3,4]
                            ¥ƒ0  | Reduce using the following as a dyad and starting with 0
                        ;  ¥     | - Concatenate the following as a dyad
                         +       |   - Add
                          Ṁ      |   - Take the maximum
                               Ṁ | Finally take the overall maximum

Python reference

Finally, here’s a TIO version of the Python reference that @pasbi included for the main problem. It reads from stdin.

Python 3 + numpy + scipy, 248 240 bytes

from scipy.optimize import*
from numpy import*
def f(b,i=0):
 for r,c in b:p=zeros(2*len(c)+1);p[-3::-2]=c;p[-1]=h=max([0,*(-fminbound(lambda x:polyval(polysub(p,d),x),0,min(s,r),full_output=1)[1]for s,d in b[:i])]);b[i][1]=p;i+=1
 return h

Try it online!

-8 bytes thanks to @xnor

The function takes a list of [radius, polynomial] pairs as input and returns the pile height.

This solution uses more or less the same algorithm as the reference code except that it does not compute the maximum using derivatives. Meanwhile, it is written using built-in numpy and scipy functions in Python. The ungolfed version is shown in the following. This serves as an alternative version of the reference code for those who wants a shorter version to capture the idea quickly.

from scipy.optimize import fminbound
import numpy as np

def compute_pile_height(bowl_data):
    for i, (r, curve) in enumerate(bowl_data):
        distances = [0]  # Initialize the distances array with 0 as the lower bound for max
        # Construct a complete polynominal coefficient array
        curve_poly = np.zeros(2 * len(curve) + 1)
        curve_poly[-3::-2] = curve
        
        # Iterate over all bowls under the current bowl
        for j in range(i):
            b_r, b_curve_poly = bowl_data[j]

            # Calculate the difference polynominal between the current bowl and bowl j
            diff = np.polysub(curve_poly, b_curve_poly)

            # Find the maximum height difference between bowl j and the current bowl in the range [0, min(b_r, r)]
            max_height_diff = -fminbound(lambda x:np.polyval(diff, x), 0, min(b_r, r), full_output=True)[1]
            distances.append(max_height_diff)

        # Compute the maximum distance as the height for the current bowl, 
        # update the polynominal using the height as the constant coefficient
        curve_poly[-1] = height = max(distances)

        # Update stored data for the current bowl
        bowl_data[i][1] = curve_poly
    return height

Try it online!

R, 451 436 bytes

function(x){x=c(x[1],x);a=rev(pmax(0,c(combn(x,2,function(y,z=sapply(y,"length<-",max(lengths(y)))){z[is.na(z)]=0
b=rep(0,2*(n=nrow(z)-1))
b[2*1:n]=e=z[-1,2]-z[-1,1]
b=b*1:(2*n)
while(!c(b,1)[1])b=b[-1]
b=rev(b)
s=`if`(length(b)>1,eigen(rbind(-(b/b[1])[-1],cbind(diag(length(b)-2),0)))$va,0)
max(outer(c(pmin(abs(s[s==abs(s)]),r<-min(z[1,])),r),2*1:n,`^`)%*%e)}))))
o={}
for(i in seq(a=x[-1])){o=c(o,max(c(0,o)+a[1:i+F]));F=F+i}
max(o)}

Try it online!

Try it online!

Broadly speaking an R port of my Jelly answer, though since base R has no function to find the roots of polynomials this is implemented using the method found in polynom::solve.polynomial.

A function taking a list of numeric vectors from top to bottom of pile.

Thanks to @RobinRyder for golfing off 15 bytes!