but those were for cartesian coordinates, and I have to work with the components x1,x2,x3 expressed in spherical ones

Whatever that means. I understand "spherical coordinates" as a set { r, theta, phi }.

The conversion to cartesian coordinates are the well known x1, x2, x3 given above.

Last not least I think the problem should be tackled in cartesian coordinates.

Consider a paraboloid "pointing" in the negative y-direction:

para = {a u Cos[v], -u^2, b u Sin[v]};

This surface is an image of a (planar) region (subset) of R2 with 0<= u < Infinity and 0<= v <= 2Pi. All points with para . para == r^2 lie on or in the surface of a sphere with radius r.

Now we can determine those u (as a function of v) which fullfill this condition, giving somewhat complicated expressions ( here given in para1) , but in principle it is quite easy to extract the spherical coordinates, at least numerically.

The follwing code gives the line of intersection (on the sphere) for different paraboloids

R = 3;
Manipulate[
 para = {a u Cos[v], -u^2, b u Sin[v]}; 
 sol1 = Flatten[Solve[para.para == R^2, u]];
 para1 = para /. # & /@ sol1;
 pl1 = ParametricPlot3D[
   para1,
   {v, 0, 2 Pi},
   PlotStyle -> {Thick, Blue}];
 pl2 = Graphics3D[Sphere[{0, 0, 0}, R]];
 Show[pl2, pl1],
 {a, -5, 5}, {b, -5, 5}]

and here we have an element of para1 and its spherical coordinates:

vec = Re[para1 /. v -> .345][[4]]
th0 = ArcCos[vec[[3]]/R]
phi0 = ArcTan[vec[[1]], vec[[2]]]
R {Cos[phi0] Sin[th0], Sin[phi0] Sin[th0], Cos[th0]}

I think

s2 = Solve[eqs2[r, \[Theta], \[Phi]] == R^2, r][[2]] // Simplify;

is not appropriate. We are looking for an intersection of the paraboloid and the sphere. All these points are elements of the spherical surface and therefore have a r-value of R, so we must not look for r.

The intersection being a line means that it is a one-dimensional entity, depending only on ONE pararmeter. So I think fixing r to the radius of the sphere we should ask for theta as function of phi or phi as function of theta and insert this in the representation of a point:

For example

s3 = Solve[eqs1[R, \[Theta], \[Phi]] == R^2, \[Phi]]
Length[s3]
plot3 = ParametricPlot3D[{x1[R, \[Theta], \[Phi]], 
    x2[R, \[Theta], \[Phi]], x3[R, \[Theta], \[Phi]]} /. 
   s3[[5, 1]], {\[Theta], 0, Pi},
  PlotStyle -> {Thick, Blue}]
Show[plot2, plot3]

or

s4 = Solve[eqs1[R, \[Theta], \[Phi]] == R^2, \[Theta]]
Length[s4]
plot4 = ParametricPlot3D[{x1[R, \[Theta], \[Phi]], 
    x2[R, \[Theta], \[Phi]], x3[R, \[Theta], \[Phi]]} /. 
   s4[[6, 1]], {\[Phi], 0, 2 Pi},
  PlotStyle -> {Thick, Blue}]
Show[plot2, plot4]

As s3 and s4 yield several solutions it seems that the intersection is given piecewise. But put together it seems that both solutions are equivalent, as it should be:

plot5= ParametricPlot3D[
{x1[R, \[Theta], \[Phi]], x2[R, \[Theta], \[Phi]], 
x3[R, \[Theta], \[Phi]]} /. # & /@ Flatten[s3],
{\[Theta], 0, Pi},
PlotStyle -> {Thick, Blue}]

plot6= ParametricPlot3D[{x1[R, \[Theta], \[Phi]], x2[R, \[Theta], \[Phi]], 
x3[R, \[Theta], \[Phi]]} /. # & /@ Flatten[Drop[s4, 3]],
{\[Phi], 0, 2 Pi}, PlotStyle -> {Thick, Blue}]

Show[plot2, plot5, plot6]

Thank you again! This solution works like a charm!

Best Regards

@Gianluca thank you so much for the answer!!!
Your solution works as intended for me, however, now I took the problem to another level and I am facing an issue that I can't solve.
In my initial problem, my second surface was a sphere, and the first one a basic paraboloid. The equation of the paraboloid was straightforward.
Now, I still have the sphere, but the paraboloid is parameterized by three numbers (in the code: e, u and v).

Hsim[r_, theta_, fi_, u_, v_] := 
  x2[r, theta, fi]^2 + u*x3[r, theta, fi]^2 + 2 v*x1[r, theta, fi];
eqs1[e_, u_, v_] := NSolve[Hsim[r, \[Theta], \[Phi], u, v] == e, r];
surface1[e_, u_, v_] := r /. eqs1[e, u, v][[2]];

I did the same thing and defined surface1 as being the solution to an equation, but this time the equation is Hsim[r,\theta,\phi,e,u,v]==e, which again I have to solve for r.

The problem is in the MeshFunction. If I try in the same manner MeshFunction -> Function[{x, y, z, r, \[Theta], \[Phi]}, Evaluate[surface1[e, up, v0p]]], I don't get the intersection curve (ellipse) between the sphere and the paraboloid.
e is just a number that I can modify in a Manipulate, to see how the intersection curve modifies with the increase if this value. The parameters u,v (denoted in the document with up, and v0p are know beforehand). What is the issue here?

PS: I even tried putting the analytic form of the surface in the MeshFunction like so:

MeshFunctions -> {Function[{x, y, z, 
    r, \[Theta], \[Phi]}, (0.3431771859697404` (7.3454135`*^7 Cos[\
\[Theta]] + 
         1.96286546`*^8 Sqrt[
          0.14003978386304503` Cos[\[Theta]]^2 + 
           4.` e Cos[\[Phi]]^2 Sin[\[Theta]]^2 + 
           0.3582649439870417` e Sin[\[Theta]]^2 \
Sin[\[Phi]]^2]))/(1.34722129`*^8 Cos[\[Phi]]^2 Sin[\[Theta]]^2 + 
       1.2066554`*^7 Sin[\[Theta]]^2 Sin[\[Phi]]^2) /. {e -> e}]}

But the lines that I get are not ok.

I attached here a simple document with the example. My first manipulate is without the MeshFunction, and that works fine. The second Manipulate is where I try to implement the intersection curves, unsuccessfully.

Attachments:

simplifiedHamilt...nb