That reminds me of Jones' formula, for which once there was a posting (in the
Application Center ? - have only filed the files ... please search, here is a copy)

# Jones Formula to compute the n-th prime number

restart:
M := proc (x::integer, y::integer)
  if x  else x-y;
  fi;
end;

JonesM := proc (n::integer) 
  local m,s,i,p,j,f,k;
  m := n*n;
  s := 1;
  for i from 1 to m do
      p:=0;
         for j from 1 to i do
             f := 1;
          for k from 1 to (j-1) do
                  f := f*k*k; 
                     f := f mod j;
              od;
           p := p+f;
         od; 
         s := s+ M(1, M(p,n));
  od;
RETURN (s);
end;

# check it
seq(JonesM(n), n=1 .. 10);

                  2, 3, 5, 7, 11, 13, 17, 19, 23, 29

Of course there is no 'practical' use for it. But it is nice :-)

Edited: use 'option remember', if you want to plot using it,
else it will recalculate at every plot point, someting like

plot('JonesM'(floor(n)), n=1 .. 10);
plot('JonesM'(floor(n)), n=1 .. 10, style=point);

Privacy would be my major concern for FB nonsense ...

If you are using Firefox, then you may look up for some add ins for privacy
improvement, mainly Ghostery (kick off sniffing), Adblock Plus (advertising),
Cookie Monster (Flash cookies) and NoScript (controlling for Java as well,
allow selectively javascript - for this site you have to allow mapleprimes
and googleapis.com because of Ajax [or even restrict to that])

All of them are free, provided through Mozilla org

And should kill your problem as a side effect, no need to argue with admins.

For me (at home) at least the following are blocked: 21 scripts and all of
http://www.ghostery.com/apps/facebook_connect,
http://www.ghostery.com/apps/facebook_social_plugins
http://www.ghostery.com/apps/omniture

May be, I should take more care for twitter ...

The last 4 of your 7 equations are linear in a1 ... a6 (and w is a way to give
a nonhomogenity for that - it will be chosen by you anyway), parametrized by t,
and t only appears in the coefficients for that linear system (rational + log
in t (so t seems to be assumed as positive).

The first 2 are 'quadratic', the 3rd is of degree 4, the coefficients are 
exp( rational + log in t ), variables are a1 ... a5 and a6 is not present.

What I would try (after convert(sys, rational) to avoid floats) is to solve the
first 4 linear equations - giving a 3 dimensional vector space, parametrized by t,
which means you have *3* degrees of freedom (besides t)

And then feed it to equation 1 - 3 (may be even to complete squares for eq 1 + 2,
i.e. a quadratic normal form and go to eq 3 in case of deeper interest).

Note however that - what way ever you go - it seems to me, that you actually have
4 variables which you are going to feed to eq 1 - 3.

Thus generically I would expect a curve as solution. It may be, that in the Real
case (which I guess you are working with) it nevertheless may result in isolated
points, but that is not clear at all.

Hope that I got it right.

Your last line means, you want to simplify subject to 1/m^2 = 0,
which probably is not what you want.

I think the appended sheet gives the proof, that in z = 1/2 the jump is +- Pi.

hirnyk_discontinuit.mws

Converting to Rationals and back to Floats works (and Maple asserts,
that this operation is exact), but will have side effects, if you also
have some symbols in the set (but that can be handeled).

  {Float(10.00)} union {Float(10.0)};
  convert(%, rational);
  evalf(%);
                                {10.}

@Markiyan Hirnyk

It is a bit ugly because of notations, conventions and handling.

  indicatorUNIT:= (x,y) -> 
    `if`( 0<=x and x <=1 and 0<=y and y<=1 and y <= -x+1,1,0);

For that it is obvious, how to handle the open 'unit triangle'.

  plot3d(indicatorUNIT, -0.2 .. 1.2, -0.2 .. 1.2, axes=boxed, 
    labels=["x","y",""], orientation=[-90,0]); #[-20,60]);

Now let be given 3 points defining a non-degenerated triangle
(not collinear is the classical notion, I think).

  B:=proc(L) # list of 3 points in the plane
  local M;
  M:=convert([L[2] - L[1], L[3]-L[1]], Matrix);
  M:=LinearAlgebra:-Transpose(M);
  return 1/M; #map(simplify, %, size);
  end proc;

  indicator:=proc(L, x,y )
  local u,v;
  B(L).Vector([x-L[1][1],y-L[1][2]]);
  u,v:= %[1], %[2];
  return indicatorUNIT(u,v);
  end proc;

Test it with L:= [[-1,-1], [1,1], [0,2]];

  check:= proc(x,y) indicator(L, x,y); end proc;

  plots:-contourplot(check, -1.5 .. 1.5, -1.5 .. 2.5, grid=[50,50], 
    filledregions=true, coloring=[white,red]);

It would be better to compute B(L) only once and it can certainly
be improved, it is a brute translation of what I sketched.

Have you looked at "save as maple input" = *.mpl?

It would also avoid the problem of classical /standard interface,
it can be read with any text editor.

L:=[seq(x(i), i=0 .. 20)]:
map('t -> type(t, integer)', L);
convert(%, set); # to see, which values occure

                                {true}

z := Complex(x, y); evalc(%);

                               x + y I

Ok, but what should diff be in a general situation, where the
function is not analytic?

g:= t -> Re(t) + Im(t)^2*I;

Then

g(z);
                                    2
                               x + y  I

and 

Physics:-diff(g(z), z);
                                  0
Physics:-diff(g(Complex(x, y)), Complex(x, y));
                                  0

Does that make sense?

# rewrites the recursion
x(n+2)=(((x(n+1)*x(n))^2+x(n)^2+x(n+1)^3))/x(n-1);
eval(%, n=n-2);

x:=proc(n)
option remember; # <--- 
(x(n-1)^2*x(n-2)^2+x(n-2)^2+x(n-1)^3)/x(n-3); 
end proc;
x(0):=x0; x(1):=x1; x(2):=x2; x(3):=x3; # initial values

Now - for example - you call it by x(8).

Have not worked it out, but I neither would use indexed variables nor arrays for
storage and also think, that 1 loop is enough, like this

for ...

  v_new = ...
  E_new = ...
  s_new = ...

  s_old = s_new, E_old = E_ new, ...

end for

Also you do not use previous values of your u, no updates while looping ?

Roughly you want

6.6e18906*n*((-2.302585093+ln((-6.*n+7.*n*m+4.*n*t+8.-10.*t-10.*m)/(-2.*n+2.*n*m+n*t-3.*m-4.*t+3.)))/(.5e6990-.1e6990*n)
 ^(7447/10000)/m)^(10000/2753)*m +
690.7755279-.3e3*ln((-6.*n+7.*n*m+4.*n*t+8.-10.*t-10.*m)/(-2.*n+2.*n*m+n*t-3.*m-4.*t+3.)) = 0.

As it was said 6.6e18906 is likely to kill numerics - where you would have to supply concrete
values for m=taum and t=taut to work at all.

And symbolically I have doubts there are usable solutions.

Best would be to think about appropriate scaling of variables, which you have to do before
arriving at your final task.

Just as thoughts

u[x](x)*(y/delta(x)-2*y^2/delta(x)^2+y^3/delta(x)^3+
   B*H*(1-3*y^2/delta(x)^2+2*y^3/delta(x)^3)/delta(x))/
  (3*B*H/delta(x)+2);
 
 Int(%^2, y=0 .. delta(x)): value(%):
 
 simplify(%);

                2                   2                         2  2
         u[x](x)  delta(x) (delta(x)  + 11 B H delta(x) + 39 B  H )
   1/105 ----------------------------------------------------------
                                               2
                           (3 B H + 2 delta(x))

find the pdfs attached maple10_sheets_pd.zip