I do not understand why you use the loop. A simple solution without any loops:

restart;

Z:=Matrix(2, symbol=z);

Z[1,1]:=1;

Z;

The command  identify does the inverse transformation:

evalf( ln(3) );

evalf( exp(5) );

identify(%%);

identify(%%);

convert(%%, ln);

I understand the question so that it is not allowed to use any packages, only the core of Maple. Solution for rows (similarly for columns):

M:=<1,  5,  6; 0,  0,  0;  2,  3, 5 >:

r:=0:

for i from 1 to op([1,1], M) do

if convert([seq(M[i,j]<>0, j=1..op([1,2],M))], `or`) then r:=r+1 fi;

od:

r;

                                           2

If I understand the problem, all the pieces are in the form of a circular sector of the same radius. The first step - add them together to get a single sector.

1) The solution of 3 sections: just cut this single sector by 3 sections into 4 equal pieces.

2) The solution of 2 sections: one section goes along the symmetry axis of the sector, the second section is perpendicular to the first one. For reasons of continuity it is always possible.

Addition. With Maple (if the radius and the angle of the sector are known) easy to find at what distance from the center of the circle to do the second cut.

 Addition 2. I downloaded the file and saw that I understood the problem incorrectly. Therefore, my solutions are the solutions of other problem. Tomleslie's solution is right.

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