surf1[a_, b_] =
ContourPlot3D[
y^2/a^2 + z^2/b^2 == 1, {x, -5, 5}, {y, -5, 5}, {z, -5, 5},
RegionFunction ->
Function[{x, y, z},
z^2/b^2 + x^2/a^2 >= 1 && z >= 0 && Abs[x] <= Min[a, b] &&
Abs[y] <= Min[a, b]], PlotPoints -> 50,
RegionBoundaryStyle -> None, Mesh -> None];
surf2[a_, b_] =
ContourPlot3D[
z^2/b^2 + x^2/a^2 == 1, {x, -5, 5}, {y, -5, 5}, {z, -5, 5},
RegionFunction ->
Function[{x, y, z},
y^2/a^2 + z^2/b^2 >= 1 && z >= 0 && Abs[x] <= Min[a, b] &&
Abs[y] <= Min[a, b]], PlotPoints -> 50,
RegionBoundaryStyle -> None, Mesh -> None];
surf[a_, b_] :=
Show[surf1[a, b], surf2[a, b], Boxed -> False, Axes -> False];
surfs = Table[surf[a, b], {a, 3, 5, .5}, {b, 3, 5, .5}];
Graphics3D[
Table[GeometricTransformation[surfs[[i, j]][[1]],
TranslationTransform[15 {i, j, 0}]], {i,
Dimensions[surfs][[1]]}, {j, Dimensions[surfs][[2]]}],
Boxed -> False, Axes -> {True, True, False}]
