restart;
with(LinearAlgebra): interface(rtablesize=infinity): with(linalg):
P:=matrix([
[1, t_1*(epsilon*t_0^2+1)*(-t_0*t_1+epsilon)*(t_1^3+epsilon*t_0)/( epsilon*t_0^2*(t_1^4-1)*(t_1^2+epsilon) ),
-t_1*(t_0^2+epsilon)*(t_0*t_1^3+1)*(t_1-t_0)/( t_0^2*(t_1^4-1)*(t_1^2+epsilon) ) ],
[1,-t_1*(t_0^2*t_1^2-1)/(t_0*(t_1^4-1)), (t_0-t_1)*(t_0*t_1^3+1)/(t_0*(t_1^4-1)) ],
[1,-( epsilon*t_0^2*t_1+t_0*(t_1^4-1)-epsilon*t_1^3 )/(t_0*(t_1^4-1)),epsilon*t_1*(t_0^2-t_1^2)/(t_0*(t_1^4-1))]
]);
E:=(t_0-t_1)*(t_0*t_1-epsilon)/(epsilon*t_0*(t_1^2+epsilon));