The equality  int(diff(f(x), x), x=a..b) = f(b) - f(a)  is not always true. For example, consider the integral  int(1/x, x=-1..2) . 

'continuous' option is needed:

Fun := Diff(f(x), x);

int(Fun, x = 2 .. 7, 'continuous');

Since your function  y->int((y⁴+by²+c)^-1/2, y)   is represented by the indefinite integral, for uniqueness, I wrote it as a definite integral  f:=y->Int((s^4+b*s^2+c)^(-1/2), s=0..y). The domain of the inverse function  y=h(x)  for the function  f(y)=x  coincides with the range of  f  when  y=0..infinity.

The procedure  h  numerically solves the problem:

restart;

h:=proc(x,b,c)

global M;

local f;

if b^2-4*c>=0 then error "Should be b^2-4*c<0" fi;

M:=Int((s^4+b*s^2+c)^(-1/2), s=0..infinity);

f:=y->Int((s^4+b*s^2+c)^(-1/2), s=0..y);

if x<0 or x>evalf(M) then error "x does not belong to the range of the function f" fi;

fsolve(f(y)=x, y=0..infinity);

end proc:

The domain of h is  the range  0 .. M .

The examples of use:

h(0.5,2,3);  

h(1.5,2,3);

evalf(M);

plot('h(x,2,3)', x=0..1, color=red, thickness=2);

The procedure  coefff  returns the coefficient of the monomial  t  in the multivariate polynomial  P (T - the list of variables of the polynomial)

coefff:=proc(P, T, t) 

local L, H, i, k:

L:=[coeffs(P, T, 'h')]: H:=[h]: k:=0:

for i from 1 to nops(H) do

if H[i]=t then k:=L[i] fi:

od:

k;

end proc:

Solution of your problem:

Id:=<3*x^2+2*x*y+x, y-x*y+3, y^2-2*x+4>:

L:=[x^2, x*y, x, y^2, y, 1]:

M:=Matrix(3,6, (i,j)->coefff(Id[i], [x,y], L[j]));

Radical of an integer is (in number theory) the product of the distinct primes which divide this integer.

The procedure for calculation:

Rad:=proc(n::integer)

ifactors(n);

`*`(seq(%[2,i,1], i=1..nops(%[2])));

end proc:

Examples:

Rad(54);

Rad(720);

     6

    30

sigma:=167:

t:=i->exp(-i^2/2/sigma^2):

S0:=l->1+2*sum(t(i), i=l..infinity):

S:=l->`if`(l>=32,evalf(S0(l)), evalf(2*add(t(i), i=l..31))+evalf(S0(32))):

S(1),  S(20),  S(31),  S(32),  S(40);  # Examples

              418.6069218, 380.6953067, 358.9442559, 356.9784188, 341.3372758

Obviously, the problem reduces to the calculation  cos(z)  and  sin(z) , if  z=1/2*arccos(1/3) . Using Euler's formula  exp(I*z)=cos(z)+I*sin(z), we obtain

expr:=expand(sin(-(1/6)*Pi+(1/2)*arccos(1/3)));

z:=1/2*arccos(1/3):

evalc(exp(I*z));  # Euler's formula in Maple

eval(expr,[cos(z)=Re(%),sin(z)=Im(%)]);

The easiest way to draw graphs of numerical solutions of differential equations - is to use  plots:-odeplot  command.

Example with 2 functions:

Sol:=dsolve({diff(y[1](x),x)=sin(y[2](x))+2*cos(y[1](x)), diff(y[2](x),x)=2*cos(y[2](x))+3*sin(y[1](x)), y[1](0)=0, y[2](0)=1}, numeric):

plots:-odeplot(Sol, [[x,y[1](x)], [x,y[2](x)]], x=0..4, color=[red,blue], thickness=2, view=[0..4, 0..4]);

Addition. If you still want to use  plot  command, you can write as in example: 

y[1]:=x->rhs(Sol(x)[2]):

plot(y[1], 0..4);

Your contour - it is an ellipse. To maintain the proportions along the axes Ox and Oy is necessary to take the same ranges. For clarity, the cutting plane  z=160   is also shown:

A := plots:-contourplot3d(4*x^2+9*y^2, x = -8.5 .. 8.5, y = -8.5 .. 8.5, contours = [160], color = red, thickness = 4, linestyle = solid):

B := plot3d(4*x^2+9*y^2, x = -8.5 .. 8.5, y = -8.5 .. 8.5, style = surface, view = -15 .. 190):

C:=plot3d(160, x=-8..8, y=-6..6, transparency=0.5):

plots:-display(A, B, C, axes = normal, orientation = [55, 55], lightmodel = light4);

Below we prove this using Maple for calculations.
First, we prove that the number  n  can not have more than 21 digits.

The equation  9^20*(n+1)=10^n  has not more than 2 roots as the intersection of a straight line with a convex curve. Find these roots and the plot (for clarity) :

Digits:=25:

fsolve(9^20*(n+1) = 10^n, n = -2 .. 0);

fsolve(9^20*(n+1) = 10^n, n = 0 .. infinity);

plot([9^20*(n+1), 10^n], n = -2 .. 22, 0..4*10^20, thickness=2);

Next, check that the number of 199 ... 9 (20 nines) satisfies the required property:

is(1+9^20*20>=10^20+10^20-1);

                           true

It remains to verify that any 21-digit number greater than 199 ... 9, does not satisfy this property. There is no need to check all the numbers up to 10^21-1. The basic idea - the sum  20th powers of the digits from the number  199 ... 9  to the number  99 ... 9  rises very slowly (just 1.05 times), so to check the property sufficiently in the small range. Changing 4 digits and doing the simple estimates we get the desired result:

restart;

evalf(9^20*21/(1+9^20*20));

evalf((1+9^20*20)/(10^20+10^20-1));

%*%%;

(2*10^20-1)*%;  # The maximum number up to which we should to check

s:=1: 

for a from 0 to 5 do

for b from 0 to 5 do

for c from 0 to 3 do

for d from 0 to 2 do

N1:=2*10^20+a*10^19+b*10^18+c*10^17+d*10^16:

N2:=2*10^20+a*10^19+b*10^18+c*10^17+d*10^16+10^16-1:

if  2^20+a^20+b^20+c^20+d^20+16*9^20 < N1 then s:=s*1 else s:=s*0 fi;

od: od: od: od:  

s;

with(numapprox):

pade(cos(x), x, [4,4]);

plot([cos(x), %], x=0..2*Pi, -1.2..1.5, color=[red,blue]);

We see that in the range  x=0..3  approximation is very good - the curves merge.

Obviously, this number cannot have 5 or more digits because 9^3*5<10000. 

Simple enumeration solves the problem:

restart;

for k to 9999 do

L:=convert(k,base,10);

s:=add(L[i]^3, i=1..nops(L));

if s>=k then n:=k fi;

od:

n;

                         1999

simplify(abs(x)-sqrt(x^2))  assuming x::real;

                                    0

This issue differs from the previous in that the right-hand side of the equation is much smaller. Therefore, the possibility appeared to estimate the  number of solutions. It turned out that solutions more than 100000000. Furthermore, among the solutions, we find the solution with minimum norm.

Finding the number of solutions:

restart;

eq:=30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

isolve(eq);

ineq:=5*_Z1+10*_Z2+25*_Z3+35*_Z4+25*_Z5+40*_Z6+5*_Z7+10*_Z8<640;

L:=coeffs(lhs(ineq));

k:=0:

for _Z1 to (640-`+`(L[2..8]))/5 do

for _Z2 to (640-5*_Z1-`+`(L[3..8]))/10 do

for _Z3 to (640-5*_Z1-10*_Z2-`+`(L[4..8]))/25 do

for _Z4 to (640-5*_Z1-10*_Z2-25*_Z3-`+`(L[5..8]))/35 do

for _Z5 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-`+`(L[6..8]))/25 do

for _Z6 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-`+`(L[7..8]))/40 do

for _Z7 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-40*_Z6-L[8])/5 do

for _Z8 to (640-5*_Z1-10*_Z2-25*_Z3-35*_Z4-25*_Z5-40*_Z6-5*_Z7)/10 do

if ineq then k:=k+1 fi;

od: od: od: od: od: od: od: od:

k;  # The  number of solutions for positive integer parameters  _Z1, _Z2, ..., _Z8

Next we find a solution with minimum norm. To do this, first we find the optimal solution in two ways without regard integer, and then varying the closest integer solutions solves the problem:

restart;

eq:=30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

F:=a^2+b^2+c^2+d^2+e^2+f^2+g^2+h^2+i^2:

Optimization[QPSolve](F,{eq}, assume=nonnegative);

extrema(F, {eq}, indets(F),'s'): evalf(convert(op(s),list));

map(t->[round(rhs(t)-4),round(rhs(t)+4)], %);

Dist:=infinity:

for a from 1 to 6 do

for b from 1 to 9 do

for c from 3 to 11 do

for d from 1 to 9 do

for e from 12 to 20 do

for f from 6 to 14 do

for g from 1 to 5 do

for h from 1 to 5 do

for i from 3 to 11 do

if eq then dist:=F;

if dist<Dist then Dist:=dist; Sol:=[a,b,c,d,e,f,g,h,i]  fi; fi;

od: od: od: od: od: od: od: od: od:

['a'=Sol[1],'b'=Sol[2],'c'=Sol[3],'d'=Sol[4],'e'=Sol[5],'f'=Sol[6],'g'=Sol[7],'h'=Sol[8],'i'=Sol[9]];  # The solution with minimal norm

Edit. Comments on the first  code was corrected. In fact, this code finds not all solutions but only part of the solutions for positive values of parameters. The total number solutions is  more than 106141495 solutions.

Example:

plot([x,x^2,x^3], x=-1..1, thickness=2, legend=[x,x^2,x^3]);

for your way:

unapply(diff(f(x), x), x);

                  x->2*x

Of course, the use of the operator  D  is more preferably.