Using this procedure, you can find the values of the unknown function at certain points and build its plot:

dy4:=diff(y(x),x):

eqn4:=dy4=sin(x*y(x)):

ic1:=y(0)=1:

ans3:=dsolve({eqn4,ic1}, y(x), numeric);

ans3(1);

plots[odeplot](ans3, [x,y(x)], x=0..10, thickness=2, numpoints=1000);

The example - all on the same plot:

A := plot3d(x^2-y^2, x = -1 .. 1, y = -1 .. 1, shading = zhue):

B := plot3d([1, y, z], y = -2 .. -3/2, z = -1 .. 1, shading = zhue):

plots[display](A, B, axes = frame, orientation = [40, 75], lightmodel=light4);

map(simplify@sum, op(sum(a*u[k]+b*v[k], k=1..n)));

combine(%);

restart;

f:=(x,y)->x*(x+y)*exp(y-x);

extrema(f(x,y), {}, {x,y}, `s`);

Points:=`s`;

Student[MultivariateCalculus][SecondDerivativeTest](f(x,y), [x,y]=[0,0]);

Student[MultivariateCalculus][SecondDerivativeTest](f(x,y), [x,y]=[1/2,-3/2]);

You have forgotten  expand  f  in the Taylor series in the neighborhood  x=4. Should be

f := x^(6*ln(x));

Digits:=15;

taylor(f, x=4, 7);

T2 := convert(%, polynom);

f_value := evalf(subs(x = 5, T2));

Let   F(x,y,t,s) = K(x,y,t,s)*h_m(x)*h_n(y)*h_p(t)*h_q(s) . 

int(F(x,y,t,s), [x=0..1, y=0..1, t=0..1, s=0..1]);

Example:

Max:=-infinity:

S:=[2, -5, 4, 0, 4, 1]:

for s in S do

if  s>Max  then  Max:=s  fi;

od:

Max;

                                                        4

To find the positive roots of your equations, I first put them to the reduced form of quadratic equation, and then find the positive roots of these equations using well-known formula, taking the plus sign in front of the square root:

eq1 := 1.6*10^(-7)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)-1.6*10^(-14)*R^2*cos(t)^2+4.2*10^(-14)*R^2-1.3+2.1*10^(-9)*R*cos(t) = 0:

eq2 := 8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.9*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t) = 0:

eq3 := 8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.2*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t) = 0:

eq4 := 2.1*10^(-9)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)+1.6*10^(-14)*R^2*cos(t)^2+2.6*10^(-14)*R^2-1.3+1.6*10^(-7)*R*cos(t) = 0:

L := [eq1, eq2, eq3, eq4]:  # List of original equations

L1 := map(collect, map(t->t/coeff(lhs(t), R^2), L), R):   # List of reduced equations

Sol := map(t->-(1/2)*coeff(lhs(t), R)+sqrt((1/4)*coeff(lhs(t), R)^2-coeff(lhs(t), R, 0)), L1):   # List of positive roots of this equations  

A := plot(Sol, t = 0 .. 2*Pi, color = [blue, red, red, green], thickness = 3, coords = polar):  #  Required curves in polar coordinates

B := plot(min(op(Sol)), t = 0 .. 2*Pi, color = yellow, coords = polar, filled = true): # The intersection of the regions

plots[display](A, B);  # All together

 Addition:   min(f(t), g(t))  -  the minimum function of  f(t)  and  g(t)

Obviously the equality  holds

x=<0.3603, 1.3279, 1.6882> = 0.3603*<1, 0, 1> + 1.3279*<0, 1, 1> 

From this equality it follows that  x  is in the rowspace  of the vectors  <1,0,1>  and  <0,1,1>

For the second and third graphics positive roots for any  t  almost are the same, so in the plot we see  not 4 but 3 graphics (blue graph covers the yellow graph):

eq1:=1.6*10^(-7)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)-1.6*10^(-14)*R^2*cos(t)^2+4.2*10^(-14)*R^2-1.3+2.1*10^(-9)*R*cos(t)=0:

eq2:=8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.9*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t)=0:

eq3:=8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.2*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t)=0:

eq4:=2.1*10^(-9)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)+1.6*10^(-14)*R^2*cos(t)^2+2.6*10^(-14)*R^2-1.3+1.6*10^(-7)*R*cos(t)=0:

L:=[eq1,eq2,eq3,eq4]:

L1:=map(collect,map(t->t/coeff(lhs(t),R^2), L),R):

plot(map(t->-coeff(lhs(t),R)/2+sqrt(coeff(lhs(t),R)^2/4-coeff(lhs(t),R,0)),L1), t=0..2*Pi, color=[red,blue,yellow,green], thickness=2);

 You can easily finish the remaining yourself.

restart;

S:=b*h^2/6: M:=10^7: fb:=M/S: fall:=2400: obj:=abs(fb-fall);

eval(obj, {b=10, h=50});  # Minimum

limit(eval(obj, {b=1/n,h=2/n}), n=infinity);  # Maximum

This works:

a:=t^2;

f:=unapply(a, t);

                                     a := t^2

                                  f := t -> t^2

Example:

f(3);

                 9

Definition of the function and its plotting for  alpha=2:

f := (x, y, alpha)->piecewise((-(1/2)*Pi <= x and x < 0 or 0 < x and x <= (1/2)*Pi) and -abs(x) <= y and y <= abs(x), (x^2+y^2)^alpha*arcsin(y/x),  x = 0 and y = 0, 0,  undefined);

plot3d(f(x, y, 2), x = -3 .. 3, y = -3 .. 3, numpoints = 25000, style = surface, axes = normal, orientation = [30, 70]);

Changing  the number  alpha  you will receive a variety of functions and their graphs.

restart;

with(plots):

B := 1.5:

u2 := r*cos(th): v2 := r*sin(th):

w1 := u1+I*v1: w2 := u2+I*v2:

z1 := evalc(w1^2): z2 := evalc(w2^2):

x1 := evalc(Re(z1)): x2 := evalc(Re(z2)):

y1 := evalc(Im(z1)): y2 := evalc(Im(z2)):

f1 := proc (theta) options operator, arrow; plot3d([-x1, -y1, v1], u1 = -6 .. 1, v1 = -B .. B, grid = [50, 50], orientation = [theta, 80], color = u1) end proc:

f2 := proc (theta) options operator, arrow; plot3d([-x2, -y2, v2], r = 1 .. 1, th = -Pi .. Pi, grid = [50, 50], orientation = [theta, 80], color = black, thickness=2) end proc:

display(seq(display(f1(10*k), f2(10*k)), k = -17 .. 18), insequence = true, axes = box, view = [-1 .. 1, -1 .. 1, -B .. B]);

eq:=diff(x(t),t,t) - 3*diff(x(t),t) + 2*x(t) = 2*t -3;

inc:=x(0)=2, D(x)(0)=4;

dsolve({eq, inc});