@Thomas Richard 

Screenshot (new session);

@rlopez 

Actually, the problem from a mathematical point of view is not complicated. It's easy to formulate it as a minimization problem.
The embarrassing fact is that good numerical methods for global minimization do not exist!

@nm 

Actually that's why integrating factors were invented. See https://en.wikipedia.org/wiki/Inexact_differential_equation

@nm The ode is supposed (implicitly) to be an expanded polynomial in y'. Your examples are of this type. If you expand the ode, it will be OK.

For the general case of obtaining

G(y',y'^2,y'^3,.....,y'^n)  = F(x,y), (actually you probably want  G(x,y',y'^2,y'^3,.....,y'^n)  = F(x,y))

the problem seems to be very complicated.

For example, (x+y(x))*(1 + diff(y(x),x)^2+x^2+1)=1  ==> (x+y(x))=1/*(1 + diff(y(x),x)^2+x^2+1)
cannot be obtained this way.

The problem reduces to write a relation F(x,y,z)=0 in an equivalent form  f(x,y)=g(x,z).
It seems to be very difficult, a solution (when it exists) is obviously not unique, and probably there is no algorithm for it. 

@gaurav_rs 

My example shows that what you describe could be the normal behavior.
This depends on the nature of your problem.

sol:=
piecewise(And(a+b = 0,a-b = 0,-m+j = 0,-j-m = 0),[[x = 0, y = y, z = 0, t = t],
[x = x, y = 0, z <> 0, t = 0], [x <> 0, y = 0, z = 0, t = 0]],And(a+b = 0,a-b =
0,-m+j = 0,-j-m <> 0),[[x <> 0, y = 0, z = 0, t = 1/2*(-j-m)/x]],And(a+b = 0,a-
b = 0,-m+j <> 0,-j-m = 0),[[x = 0, y <> 0, z = 1/2*(-m+j)/y, t = 0]],And(a+b =
0,a-b <> 0,-m+j = 0,-j-m = 0),[[x = 0, y = 0, z <> 0, t = 1/2*(a-b)/z]],And(a+b
= 0,a-b <> 0,-m+j = 0,-j-m <> 0),[[x <> 0, y = 0, z = -x*(a-b)/(j+m), t = 1/2*(
-j-m)/x]],And(a+b = 0,a-b <> 0,-m+j <> 0,-j-m = 0),[[x = 0, y <> 0, z = 1/2*(-m
+j)/y, t = (a-b)*y/(-m+j)]],Or(And(a+b <> 0,a-b = 0,-m+j = 0,-j-m = 0),And(a+b
<> 0,a-b = 0,-m+j = 0,-j-m <> 0)),[[x <> 0, y = 1/2*(a+b)/x, z = 0, t = 1/2*(-j
-m)/x]],Or(And(a+b <> 0,a-b = 0,-m+j <> 0,-j-m = 0),And(a+b <> 0,a-b <> 0,-m+j
<> 0,-j-m = (a-b)*(a+b)/(-m+j))),[[x <> 0, y = 1/2*(a+b)/x, z = (-m+j)*x/(a+b),
t = 1/2/(-m+j)/x*(a^2-b^2)]],[]);

I'll try later (probably) to post a method to obtain this.

Edited.

@Kitonum 

You have a missing "-" in Sys. As I have mentioned, with u^(2/3) the conclusion is not true.

@Carl Love 

Unfortunately eliminate does not work properly here.
For example, a=b is not enough to have solutions.
So, the "general"  solution {...}  (z<>0)  is valid for   a^2 <> b^2  and  j^2 <> m^2  and  a^2-b^2+j^2-m^2 = 0  

Then the cases  a=b etc  must be considered separately.

A more reliable command is
SolveTools:-SemiAlgebraic(sys,[x,y,z,t]);   # (*)
but it is very slow in this case and I had not enough patience.
But the simpler system

SolveTools:-SemiAlgebraic([t*x=a, y*z=b, t*z=c, x*y=d], [x,y,z,t]);
works, and seems to find all the cases (for reals, it considers >0, <0, =0 ...).
Note that actually the original system can be reduced to (*).

You recently posted a similar question about the Galerkin finite element method.
After you received an answer and a few comments, you deleted it.
A potential respondent should be warned.

@mmcdara 

I don't know what you are talking about. I have not modified anything.
What insult? By without detours I ment a standalone formula, not depending on some randomly generated objects.
If you don't like/want my answers, just say it and I won't bother you in the future.

@waseem 

The method in the article works for linear problems. For nonlinear ones there are lots of complications.
 

@waseem 

Why would you want this for such a simple problem which is solved by Maple at once?
The authors had a strong motivation: they needed a published paper.

@mmcdara 

I find the Abramowitz-Stegun's 26.2.23   correct:  |eps(p)| < 4.5*10^(-4)  for 0<p<=1/2.
You should simply give your alternative better form (without detours) in order to compare.
Note that A-S formula obtains uniform best approx (with Remez). You obtain some L^2 approximation; it cannot be better in the worst case ( ||u||_2 <= ||u||_infinity ).

@Rouben Rostamian  
Yes, it works, not very intuitive (without a plot) though.

restart;
a:=1/3: h:=1/20: d:=1/4:
M:= # Moebius
<(1+u*a * cos(t/2))*cos(t), (1+u*a * cos(t/2))*sin(t), a*u * sin(t/2)>:
M_u, M_t := <diff(M,u)>, <diff(M,t)>:
NN:=LinearAlgebra:-`&x`(M_u,M_t):
N:=LinearAlgebra:-Normalize(NN,2):
p1:=plot3d( M+h*N, t=0..2*Pi, u=-1..1, scaling=constrained, color=orange):
p2:=plot3d( M-h*N, t=0..2*Pi, u=-1..1, scaling=constrained, color=green):
narr:=10:
q1:=seq( plots:-arrow(eval(M+h*N, [u=0,t=k*Pi/narr]),eval(d*N, [u=0,t=k*Pi/narr]), width=1/20, color=red), k=0..2*narr-1):
q2:=seq( plots:-arrow(eval(M-h*N, [u=0,t=k*Pi/narr]),-eval(d*N, [u=0,t=k*Pi/narr]), width=1/20, color=blue), k=0..2*narr-1):
plots:-display(p1,p2,q1,q2);

My intuition refuses to work about the normals for the Mobius strip.
Mobius strip is not orientable; the normal vector is not continuous. I cannot imagine the "correct" exterior normals for the two "translated" strips.  Isn't the 3D printer going to be confused?