You can use the Netwon-Rapshon method that is very easy to implement.
I coded it with Scilab:
Code:
function z=H(x,y), z=tan(x)*tan(x*y/(1-y))-5, endfunction // Function to be solved
eps=1e-4; // Wanted accuracy
i=1;
for yval=0.01:0.01:0.99, // y variable SPAN
x0=0.1; // first guess
while abs(H(x0,yval))>eps do // Iterate until the solution is inside the accuracy
x0 = x0 - H(x0,yval)/derivative(list(H,yval),x0);
end
solution(i,1) = x0; // Solution x
solution(i,2) = yval; // relative to this value of y
i = i+1; // increases index
end
try to use Wolfram Mathematica
for example look - for y(x)=sin(x):
http://www.wolframalpha.com/input/?i=solve%28+tan%28x%29+tan%28%28sin%28x%29%2F%281-sin%28x%29%29%29+x%29-5+%3D+0%2C+x%2C{0%2C1}%29