Creating and using an explicit piecewise function in a convenient way

I have a set of data points that define a function in the form

curvePts = {data1,data2,...}


where

data1 = {{x11,y11},{x12,y12},...}
data2 = {{x21,y21},{x22,y22},...}
...


For example

pts = {Table[{x, Tan[x]}, {x, 0, \[Pi]/2, 0.01}],
Table[{x, Tan[x]}, {x, 2 \[Pi], 5 \[Pi]/2, 0.01}]}


Now, I want to construct a piecewise function that interpolates pts, but I need an explicit representation.

In that respect, I've built the following function:

pPieces[pts_, n_] := Function[\[FormalX], Module[{noCurves, resN, tmpParts, polis,
domains},
noCurves = Length@pts;
resN = Table[Mod[Length@pts[[i]], n + 1], {i, 1, noCurves}];
tmpParts = Flatten[Table[Partition[pts[[i]], n + 1, n, {1, 2}, {}], {i, 1,
noCurves}], 1];
polis = Table[Fit[tmpParts[[i]], Table[y^k, {k, 0, Length@tmpParts[[i]] - 1}], y],
{i, 1, Length@tmpParts}];
domains = tmpParts[[All, {1, -1}]][[All, {1, -1}, 1]];
Piecewise[MapThread[{#1, #2[[1]] <= y < #2[[2]]} &, {polis, domains}],
Indeterminate] /. y -> \[FormalX]
]
];


which, in a nutshell, determines the number of disjoint data, then partitions the data into $n+1$ sets of points to do a $n$-th orded polynomial interpolation (leaving the end points with the appropriate $m$-th order polynomial approximation $(m<n)$), and joins it by using PieceWise. (Why I'm using \[FormalX] and replacement rules is more because of ignorance than knowledge).

My ultimate goal is to use NSolve with this explicit piecewise polynomial function. This in fact works with my naive example:

AbsoluteTiming[NSolve[pPieces[pts, 3]@x - 10 == 0, x, Reals]]

(* {0.098064, {{x -> 1.47112}, {x -> 7.7543}}} *)


Of course, there is no need for me to do all this with the example I've provided; in reality I'm working with very sinuous data with lots of points and I'm trying to solve a 37 x 37 system of (transcendental) equations, and I need all the real solutions (the main reason behind all the problem). If I solve a simple 2 x 2 system where pts1 = 885 and pts2 = 1027, it takes forever.

Here is an example of the real data, and a simple 2x2 system (computed in an i7):

AbsoluteTiming[NSolve[{215.879 x - 182.66 y + D[pPieces[pts1, 20]@x, x] ==
0, -163.841 x + 139.778 y + D[pPieces[pts2, 20]@y, y] == 0}, {x, y}, Reals]]

(* {1888.208001, {{x -> -3.46927, y -> -4.06635}, {x -> -0.675128, y -> -0.79221},
{x -> -0.141648, y -> -0.165793}, {x -> 0.032764, y -> 0.0397147},
{x -> 0.405087, y -> 0.4732}}} *)


where n = 20 has been chosen arbitrarily (in order to speed up computation, allegedly).

A note: My full problem requires me to find all real solutions for a large system of $k \times k$ coupled transcendental equations of the form

$$\vec{p_i} \cdot \vec{z} + f_i(z_i) = 0, \quad 1 \le i \le k, \quad \vec{z} \in \Omega \subset \mathbb{R}^k,$$

where $\vec{p_i}$ are constant vectors, and the $f$'s are real functions extracted from datasets of the type pts1 and pts2, evaluated at the $i$-th component of $\vec{z}$. I know, and I've been using FindRoot, Interpolation and the custom FindAllCrossings2D and FindAllCrossings3D (refs here and here) to attack parts of the problem, but I haven't been able to work with the full system, and this is my $n$-th approach to solve it. Any suggestions on how to tackle it will be very appreciated.

EDIT I've calculated the time to run the code with respect to $n$, as can be seen in the figure ($n=4$ being optimal). It's clear that a compromise must be found between the number of points taken to do the Fit, due the complexity of the domain in NSolve. I'm guessing some sort of adaptive algorithm to produce different fits for parts of the data, in order to maintain fidelity, but also to reduce the number of parts for the PieceWise function must be developed. Any ideas in this regard?

$\hskip1.7in$

-