@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?

@mmcdara 

It works in Maple 2019. But not properly for Digits>15 (try 25 and 20). Probably Carl will fix this.

Unfortunately you have not included a corresponding Your_NormInv to compare the accuracy.

@mmcdara 

It works in newer versions.
For Maple 2015 just replace z1=z  with  abs(z1-z)<10^(-Digits+1)   (or something similar).

@DJJerome1976 

AFAIK the Risch algorithm is not fully implemented in Maple.
You can see the methods which are tried using
infolevel[int]:=5;

BTW, rewriting just a bit the integrand, Maple integrates easily:

int(sin(x)^(1/3)*(1-sin(x)^2)*cos(x), x);
   

@nm 

Take

ode:=y(x)*diff(y(x),x) - y(x) = 0:
dsolve(ode);

   y(x) = 0, y(x) = x + _C1

This is obviously the correct answer (two solutions; actually there are other solutions by combining these on intervals).

If we solve ode wrt y'  ==>  diff(y(x),x) = 1
[obtained  dividing by y (supposed to be <> 0; that is how solve works)].

Of course, dsolve( diff(y(x),x) = 1 ) ==>  y(x) = x + _C1,
so, y=0 disappears. But y(x)=0 satisfies (of course) the ode.