Maximize objective function in intersection of ellipses

ellipses.py

import numpy as np
from matplotlib import pyplot as plt

nelips = 3
ndim = 2
all_axis = []
for i in range(nelips):
    matrix = np.random.rand(ndim, ndim)
    matrix += np.transpose(matrix)
    eigvals, eigvecs = np.linalg.eigh(matrix)
    lenghts = np.random.uniform(1, 5, ndim)
    axis = np.array([li * vi for li, vi in zip(lenghts, eigvecs.T)])
    all_axis.append(axis)

matrices = []
for i, axis in enumerate(all_axis):
    diag = np.diag(np.dot(axis, np.transpose(axis)))
    matrix = np.array([li/di for di, li in zip(diag, axis)])
    matrix = np.dot(np.transpose(matrix), matrix)
    matrices.append(matrix)
    print(f"axis {i} = ")
    print(axis)
    print(f"matrix {i} = ")
    print(matrix)

thetas = np.linspace(0, 2*np.pi, 129)
for i, axis in enumerate(all_axis):
    points = np.tensordot(np.cos(thetas), axis[0], axes=0)
    points += np.tensordot(np.sin(thetas), axis[1], axes=0)
    xvals, yvals = np.transpose(points)
    plt.plot(xvals, yvals)

plt.show()