That should be
fprintf(output, "%f\t %f\t \n", x, evalf(b)) ;

Thumb if You like.

There is an entire worksheet in that one line. It is recoverable I believe.

I've selected the entire contents of the sheet and converted it to Plain Text.
Result attached. corrupted-5.mw

Thumb if You like.

This curve is well known, parameterizable curve. You can do the following:

VR := -3 .. 3, -3 .. 3;

-3*a*x*y + x^3 + y^3 = 0;
Q1 := algcurves:-parametrization((lhs - rhs)(%), x, y, t);
plot(eval([op(Q1), t = 0 .. infinity], a = 1));
plot(eval([op(Q1), t = -1 .. 0], a = 1));
plot(eval([op(Q1), t = -infinity .. -1], a = 1));
plot(eval([seq([op(Q1), t = R], R = [-infinity .. -1, -1 .. 0, 0 .. infinity])], a = 1), color = red, view = [VR]);

x + y + a = 0;
Q2 := algcurves:-parametrization((lhs - rhs)(%), x, y, t);
plot(eval([op(Q2), t = -infinity .. infinity], a = 1), color = blue, view = [VR]);

plots[display](
plot(eval([seq([op(Q1), t = R], R = [-infinity .. -1, -1 .. 0, 0 .. infinity])], a = 1), color = red, view = [VR]), 
plot(eval([op(Q2), t = -infinity .. infinity], a = 1), color = blue, linestyle=dash, view = [VR])
);

See if this helps:
Go to:
Tools > Options > Display
and checkmark : Always insert new execution group after executing

Inside of f1 you could use:
if is(x, EvalfableProp) = true then "Your Test Goes Here" else 'procname(args)' fi
also define
D(f1):=0;

If you can anticipate the answer then this may be fastest method:
mtaylor(f1, [x, y], 24);
#or poisson(f1, [y], 24);
gfun[seriestolist](series(%, y, 24)):
gfun[guessgf](%, y, [ogf]);

a:=map(`^`, [8, 9 ,9 ,7 ,9 ,10 ,5], -1) mod 11;

Using Your suggestion I managed the code below. Nothing else worked.

G := (x, X) -> map((f, x, X) -> f(x) = subs(exp(I*x) = X, expand(convert(f(x),exp))), [sin, cos], x, X);
H := (x, X) -> X = convert(exp(I*x), trig);
expand(evala(subs(G(x, X), G(y, Y), ((expand@numer)/(expand@denom))(f1))));
simplify(subs(H(x, X), H(y, Y), %), trig);

Your desired solution differs from the general one as illustrated:
map[3](applyop, S->remove(evalb, S), 1, [eliminate(eqns union {theta[2, 1] = 1, theta[3, 0] = beta}, fc)]);

It's only a subcase (or so I believe).

Thumb if You like.

alias(mt=LinearAlgebra[Transpose], mht=LinearAlgebra[HermitianTranspose]);
map(radnormal, evalm(mht(F) &* F));
map(radnormal, evalm(mt(F) &* F));

Is this the issue:

alias(mt = LinearAlgebra[Transpose], mht = LinearAlgebra[HermitianTranspose]);
map(radnormal, evalm(mht(F) &* F));
map(radnormal, evalm(mt(F) &* F));

Thumb if You like.

Maybe like this:

Edit: added one missing * for ky4

 RK2skritt := proc(FR::procedure, xo, yo, vxo, vyo, h)
local x, y, vx, vy, r, kx1, kx2, kx3, kx4, ky1, ky2, ky3, ky4, lx1, lx2, lx3, 
  lx4, ly1, ly2, ly3, ly4, tmp;
  r := sqrt(xo*xo + yo*yo);
  tmp := - (h*FR(r))/(r);
  lx1 := h*vxo;
  ly1 := h*vyo;
  kx1 := tmp*xo;
  ky1 := tmp*yo;
  lx2 := h*(vxo + 0.5*kx1);
  ly2 := h*(vyo + 0.5*ky1);
  r := sqrt((xo + 0.5*lx1)^2 + (yo + 0.5*ly1)^2);
  tmp := - (h*FR(r))/(r);
  kx2 := tmp*(xo + 0.5*lx1);
  ky2 := tmp*(yo + 0.5*ly1);
  lx3 := h*(vxo + 0.5*kx2);
  ly3 := h*(vyo + 0.5*ky2);
  r := sqrt((xo + 0.5*lx2)^2 + (yo + 0.5*ly2)^2);
  tmp := - (h*FR(r))/(r);
  kx3 := tmp*(xo + 0.5*lx2);
  ky3 := tmp*(yo + 0.5*ly2);
  lx4 := h*(vxo + kx3);
  ly4 := h*(vyo + ky3);
  r := sqrt((xo + lx3)^2 + (yo + ly3)^2);
  tmp := - (h*FR(r))/(r);
  kx4 := tmp*(xo + lx3);
  ky4 := tmp*(yo + ly3);
  x := xo + 1/6*lx1 + 1/3*lx2 + 1/3*lx3 + 1/6*lx4;
  y := yo + 1/6*ly1 + 1/3*ly2 + 1/3*ly3 + 1/6*ly4;
  vx := vxo + 1/6*kx1 + 1/3*kx2 + 1/3*kx3 + 1/6*kx4;
  vy := vyo + 1/6*ky1 + 1/3*ky2 + 1/3*ky3 + 1/6*ky4;
  [x, y, vx, vy];
end proc;

Thumb if You like.

Th ecode in text format:
fsolve(3.*10^6 = sqrt((-4.2*10^14*(-1.*10^7*u*cos(u) + 3.84*10^8)/((-1.*10^7*u*cos(u) + 3.84*10^8)^2 + 1.*10^14*u^2*sin(u)^2)^(3/2) + 4.2*10^21*u*cos(u)/((-1.*10^7*u*cos(u) + 3.84*10^8)^2 + 1.*10^14*u^2*sin(u)^2)^(3/2))^2 + 7.056*10^43*u^2*sin(u)^2/((-1.*10^7*u*cos(u) + 3.84*10^8)^2 + 1.*10^14*u^2*sin(u)^2)^3), u, u = 0 .. 2*Pi);

Have You tried:
map(assign, Sol);

You must be observant of valuation rules, so:

for i from 1 by 1 while i <= nops(Sol) do
proc(x,y) assign(x, y) end(op(op(i, Sol)))
end do;

Thumb if You like.

Lookup ?igcdex.
If you need the code then you can read how it was coded with: showstat(igcdex);.

Thumb if You like.