Parameterization.mw
Maple15.
Numerical parameterization with an auxiliary surface. Always works in the local case and sometimes “global”.
Obtain the curve of intersection of the two surfaces, and then move it in any direction along the surface. Numerate and remember the point on the curve in each of its position. You can get a MxN matrix and use it as the index parameters of the surface.
Examples in the form of animation.
https://vk.com/doc242471809_375450085
https://vk.com/doc242471809_375799361
https://vk.com/doc242471809_375837573
https://vk.com/doc242471809_375682684  

   I think you would be able to help yourself consider "a" as a function of a = a (x, y). Approximate considerations in text.

Fa.mw    

Yuri Nikolaevich, using  print() suit you?
PRICONV.mw

     Numerically for the curve given any kind of equations. With geom3d  you construct the plane equation for any three points of the curve that does not lie on one line. In a cycle of successively take point of the curve, substitute them in the equation of the plane and monitor the absolute value of the discrepancy.

Many choices on the same basis.
PAR.mw

restart:
a := [x+1, x+2, x+3, x+4];
a := convert(a, set);
a :=minus(a, {x+2});
a := convert(a, list);

restart;
a := [x+1, x+2, x+3, x+4];
a := subs(x+2 = NULL, a);
nops(a);
op(2, a);

restart: with(plots):
t := [1, 2, 3];
f(t[1]) := [3, 4]; f(t[2]) := [11, 12]; f(t[3]) := [41, 1];
pointplot([f(t[1]), f(t[2]), f(t[3])], color = RGB(7, .3, 4), style = line, symbol = solidcircle, thickness = 5);
pointplot3d([1, op(f(t[1])), 2, op(f(t[2])), 3, op(f(t[3]))], color = RGB(7, .3, 4), style = line, symbol = solidcircle, thickness = 5);

combine(-ln(x)+ln(y), symbolic);

restart;
f := x1^2/(x2^3*x3^2);
op(1, op(2, f))^sign(op(2, op(2, f)));
op(1, op(3, f))^sign(op(2, op(3, f)));
f := algsubs(1/x2 = x2b, f);
f := algsubs(1/x3 = x3b, f);

restart;

(diff(f(x), x))/(diff(ln(x), x));

for example:

restart;

f := sin(x);

(diff(f, x))/(diff(ln(x), x));

Skeptik18(_for_d1.mw

 For a start point “a” = 0.5 any number of solutions. It depends on the "smax".

   1, [(0.8013209420000008)], 2.70183810879842667*10^-8                                                            

   2, [(1.0938038328000006)], 1.75716205141895898*10^-7                                                        

      3, [(1.5165511908)], 3.19181676505797540*10^-8                                                        

   4, [(1.9061358998000002)], 4.58627091859398206*10^-9                                                           

   5, [(2.1833650214000007)], 4.36816931514982798*10^-8                                                            

   6, [(2.5877177300000005)], 1.71876049503971729*10^-7                                                            

    7, [(2.937538889999998)], 9.29630750157173224*10^-9                                                           

    8, [(3.219756611999996)], 2.01128799837135830*10^-7                                                            

   9, [(3.6180527559999964)], 4.55144026911824540*10^-8                                                           

   10, [(3.953152181999997)], 3.39347172584325562*10^-8                                                            

   11, [(4.239272189999994)], 2.57192325658905930*10^-7                                                           

   12, [(4.635025125999992)], 3.01111789280383846*10^-7                                                         

   13, [(4.962516473999992)], 9.20600748965938465 10^-8                                                            

   14, [(5.2514410039999975)], 2.19850920579744980*10^-7                                                           

   15, [(5.645899321999999)], 1.82685308214303177*10^-7                                                         

   16, [5.968760335999992)], 2.30153605063065925*10^-8                                                          

   17, [6.2597571959999945)], 2.33606770816408017*10^-7                                                           

   18, [(6.653469225999989)], 1.16942534766906194*10^-7                                                           

    19, [(6.97322109999999)], 7.71696275769784279*10^-8                                                          

   20, [(7.265801703999997)], 4.86139439814792240*10^-8                                                           

   21, [(7.659044969999998)], 4.08131656026711200*10^-7                                                          

   22, [(7.976567196000005)], 6.48611420128730742*10^-8                                                    

   23, [(8.270394148000008)], 1.66075897922723926*10^-7                                                           

   24, [(8.663323925999999)], 2.60858493916771295*10^-7                                                      

      25, [(8.979170028)], 7.82206077687419565*10^-8                                                           

   26, [(9.274002000000005)], 1.76298078358172460*10^-8                                                           

   27, [(9.666711878000001)], 5.96099682503847818*10^-7                                                        

      28,[(9.981252484)], 2.84982932807764656*10^-8                                                            

   29, [(10.276911416000011)], 2.49167541710448860*10^-7 

...

The Draghilev method. Read, for example: http://www.mapleprimes.com/posts/145360-The-Dragilev-Method-1-Some-Mathematical

restart;
Digits := 30; 2^29.403243784;

For example:

restart:
nn := nextprime(10^100); zz := 1;
for ii from 0 to 100000 do zz := `mod`(zz^2+1, nn); if `or`(ii > 99997, zz = 66388502) then print("ii=", ii, "zz=", zz) end if end do:

Blue - the denominator sin(x+(1/3)*Pi-theta) = 0.
theta, I think, has a period of Pi, and x has a period of Pi / 3. The solution is obtained by Draghilev method. This numerical solution of ordinary differential equations with initial conditions theta (0) = 0, x (0) = Pi / 3.

METHOD(n-1)2d.mw