The addition and subtraction of the same multiple of x in the last equation (and other absence of automatic simplification) is indication of a cut&paste job.

It remains to be seen whether the OP might have made a transcription mistake, or might want a best fit (LeastSquares, say, easy enough after using GenerateMatrix).

@Gillee On my 64bit Linux the Threads:-Map is somewhat faster than the map instances, in Maple 2018.2 and 2020.1.

Keep in mind that n=100 is small, and Threads:-Map often needs a large number of items before the thread management overhead is small compared to the overall computation time. With n=100 there might not be any benefit, and the relative performance situation could even be switched on some platforms. That's not unreasonable.

You could also run with n=1000 or n=10000 and iterations=1. For n=10000 I see Threads:-Map being about twice as fast as map.

I should probably also have used a larger n to illustrate. The Original Poster had size equivalent to n=100 in his Question.

How accurately do you want to know evalf(R5)?

@radaar I wrote the FFG procedure so that I could Compile it. The Compiler will translate an `add` call to a loop anyway, so I doubt it'd improve much when compiled.

Moreover,  I've seen problems where the compiled translation messes up (casting floats to integers), so that's another reason I wrote out the loop myself.

What do you mean by, "I found out that my exact results are not correct"?

ps. I deleted one duplicate posting of this question. If you want to include a worksheet then attach it to a Comment here, instead of unhelpfully posting it in a new and separate Question thread.

You can check the exact result symbolically:

restart;
de:=diff(y(x), x) = 2*(1+x)-y(x):
Y:=dsolve({de, y(2) = 5});
                                                         
                          exp(-x)
        Y := y(x) = 2 x + -------
                          exp(-2)

eval(Y,x=2);

                y(2) = 5

evalb(eval(de, Y));

                  true

Are you trying to say that you had difficulty with your numeric computations?

@sepbergliaffa I didn't state that you ought to use just simplify(...,radical). I was trying to point out the difference between it as simplify(...,radical,symbolic), which is germane here.

The use of assumptions, and aspects of the symbolic option, is mentioned in the Help. That's why I referred to the Help.

@panke The code has to be such that the plots would be printed, or else it won't get sent to either the worksheet or any file.

I mentioned before (in two responses!) that you might want to wrap with a print call (if the plots wouldn't normally be printed due to further nesting, being inside a procedure, etc...).

If you terminate the loop with a full colon (or nest further, etc) then the plots wouldn't get printed inline either.That's why I suggested print(plot(...)) which you ignored when you used display(...).

@aoakindele What do you intend by  (D@@(F))(0)  in the boundary conditions?

@aoakindele 

You have a term
   F*(diff(theta(eta), eta))
in your DEs. Did you intend it to be, say,
   F(eta)*(diff(theta(eta), eta))

The error message states that "F(eta) and F cannot both appear in the given ODE". Doesn't that explain the problem?

@Carl Love He has inf=1 in the upper BCs (which happens to exceed the plotting range), but it's not clear yet whether he'll eventually want a greater value and run into convergence issues.

@panke Why don't you upload and attach your worksheet?

(Btw, why not utilize print around display?) [edit] Why ignore my earlier advice?

@pik1432 If you had asked that question about structure and representation on its own then I might just point you to the Help pages on:
  dismantle
  ToInert
and the Programming Guide sections:
  Internal Representation
  Types and Operands

But given your original Question in this thread, I'll explain some of the steps, from an "operand" point of view.

Download op_example.mw

@Katatonia I suggest that you upload and attach your complete code and worksheet to a Reply.

(The failure of the numeric integration involving Y5 may be just a matter of working precision and tolerance, but why guess?)

@pik1432 One of my points is that it's easy to simply type this in:

M__a*I__a*combine(sin(omega*t+alpha)*sin(omega*t+phi));

       M__a I__a (1/2 cos(alpha - phi) - 1/2 cos(2 omega t + alpha + phi))

But doing so (and the same for similar approaches where you have to manually separate or introduce the M__a*I__a factor) is rather ad hoc.

I'd like an ARM port of the kernel and command-line interface.