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.

@Rouben Rostamian  

OP wants an asymptotic expansion for G0(s). He should use asympt(G0(s),s)

The question does not seem to have much sense.

@9009134 

# K[r, s]
for p from 1 to 3 do  
for  s from 1 to 3  do
ff := L[p, s](r, theta, phi)*F[p, 1](r, theta):
Aij:=[indets([F[p,s](r,theta),ff],indexed)[]];
# print(indets=Aij);
g1:=expand(ff, Aij, distributed);
g2:=coeffs(g1, Aij, 'T'):
C1:=int~([g2], theta = 0 .. 2*Pi, r = 0.5 .. 1, epsilon=1e-8, numeric):
kk:=add(C1 .~ T);
kkk := evalf( int(F[p, 1](r, theta)^2, theta = 0 .. 2*Pi, r = 0.5 .. 1) );
k[p, s] := kk/kkk;
print([p, s] = %);
od
 od;

@Carl Love 

Unfortunately, as I see, sw_zgeevx_  also fails sometimes.

The bug seems to occur only for multiple (repeated) eigenvalues in the non-symmetric case, so, in practice should be extremely rare.
The examples were found using e.g.

LinearAlgebra:-CompanionMatrix((x-15)^3, x)^+;