gh-35826: Fixes to support numpy 1.25.0
    
### 📚 Description

Mostly due to a deprecation warning:
`Conversion of an array with ndim > 0 to a scalar is deprecated, and
will error in future. Ensure you extract a single element from your
array before performing this operation. (Deprecated NumPy 1.25.)`.

Also fix:
 - a change of output in `numpy.linalg.eig()`
 - a change in an exception class

See: #35081

<!-- Please provide a concise, informative and self-explanatory title.
-->
<!-- Don't put issue numbers in the title. Put it in the Description
below. -->
<!-- For example, instead of "Fixes #12345", use "Add a new method to
multiply two integers" -->


<!-- Describe your changes here in detail. -->
<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. It should be `[x]` not `[x
]`. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.


<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #35826
Reported by: Gonzalo Tornaría
Reviewer(s):

Author: Release Manager

src/sage/calculus/desolvers.py

@@ -1598,7 +1598,7 @@ def desolve_odeint(des, ics, times, dvars, ivar=None, compute_jac=False, args=()
         sage: ic=epsilon
         sage: t=srange(0,2/epsilon,1)
         sage: sol=desolve_odeint(f,ic,t,y,rtol=1e-9,atol=1e-10,compute_jac=True)
-        sage: p=points(zip(t,sol))
+        sage: p=points(zip(t,sol[:,0]))
         sage: p.show()
 
     Another stiff system with some optional parameters with no
@@ -1637,7 +1637,7 @@ def desolve_odeint_inner(ivar):
                 J = fast_float(J, dvar, ivar)
 
                 def Dfun(y, t):
-                    return [J(y, t)]
+                    return [J(y.item(), t)]
 
         # n-dimensional systems:
         else:

src/sage/matrix/matrix2.pyx

@@ -430,12 +430,12 @@ cdef class Matrix(Matrix1):
             try:
                 return self.transpose().solve_right(B, check=check)
             except ValueError as e:
-                raise ValueError(str(e).replace('row', 'column'))
+                raise e.__class__(str(e).replace('row', 'column'))
         else:
             try:
                 return self.transpose().solve_right(B.transpose(), check=check).transpose()
             except ValueError as e:
-                raise ValueError(str(e).replace('row', 'column'))
+                raise e.__class__(str(e).replace('row', 'column'))
 
     def solve_right(self, B, check=True):
         r"""

src/sage/matrix/matrix_numpy_dense.pyx

@@ -382,8 +382,9 @@ cdef class Matrix_numpy_dense(Matrix_dense):
             sage: m = matrix(RDF,[[1,2],[3,4]])
             sage: n = m.numpy()
             sage: import numpy
-            sage: numpy.linalg.eig(n)
-            (array([-0.37228132,  5.37228132]), array([[-0.82456484, -0.41597356],
+            sage: tuple(numpy.linalg.eig(n))
+            (array([-0.37228132,  5.37228132]),
+             array([[-0.82456484, -0.41597356],
                    [ 0.56576746, -0.90937671]]))
             sage: m = matrix(RDF, 2, range(6)); m
             [0.0 1.0 2.0]

src/sage/plot/plot3d/list_plot3d.py

@@ -602,7 +602,7 @@ def g(x, y):
         from .parametric_surface import ParametricSurface
 
         def g(x, y):
-            z = f([x, y])
+            z = f([x, y]).item()
             return (x, y, z)
         G = ParametricSurface(g, (list(numpy.r_[xmin:xmax:num_points * j]),
                                   list(numpy.r_[ymin:ymax:num_points * j])),

src/sage/plot/plot3d/plot3d.py

@@ -378,7 +378,7 @@ def to_cartesian(self, func, params=None):
             ....: [ 0.16763356,  0.19993708,  0.31403568,  0.47359696, 0.55282422],
             ....: [ 0.16763356,  0.25683223,  0.16649297,  0.10594339, 0.55282422]])
             sage: import scipy.interpolate
-            sage: f=scipy.interpolate.RectBivariateSpline(v_phi,v_theta,m_r)
+            sage: f=scipy.interpolate.RectBivariateSpline(v_phi,v_theta,m_r).ev
             sage: spherical_plot3d(f,(0,2*pi),(0,pi))
             Graphics3d Object