First, I believe that you mean to assume that $T$ is affine on the intervals $I_1$ and $I_2$ for, otherwise, Moran's formula need not hold. Assuming so, your question essentially boils down to how we can deal with the function $s(x,y)$ numerically. Given $x$ and $y$, the equation
$$x^s + y^s = 1$$
can easily be solved for $s$ using Newton's method and there are plenty of software tools that can be used to do so. Here is a simple way to define $s(x,y)$ and plot it using Mathematica.
sol[x_, y_] := s /. First[
FindRoot[x^s + y^s == 1, {s, 1}]];
Plot3D[sol[x, y], {x, 0, 1}, {y, 0, 1},
RegionFunction -> Function[{x, y}, x + y < 1]]
If we're a little more clever and force the graph to have the correct values on the boundary, we get something like the following:
Once we have a reliable way to compute $s(x,y)$, there are several options to compute its second derivatives. There are difference quotients that approximate the second derivatives for example. Alternatively, we might approximate $s(x,y)$ with a piecewise smooth interpolation and then compute the second derivatives based on this interpolation. Taking this second approach, I was able to construct the following graphs of $s_{xx}$ and $s_{xy}$. Of course, $s_{yx}=s_{xy}$ and $s_{yy}(x,y)=s_{xx}(y,x)$ by symmetry.
There's plenty more that could be done but it's hard to say exactly what more should be done without more specifics.