f1 can be integrated formally, but it seems f2 cannot.
(the assumptions I use come from the kc and kh ranges you give for f1).

Download f1_and_f2.mw

The simple command

res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars})):

works perfectly well and your code produces something which seems correct.
Did you think that  tolerance is an option of eval? In this cas the answer is NO: eval has only 1 or 2 arguments (look to the corresponding help page)

Had it been the case, the lines

tolerance := 1e-6:
res := eval(<vars>, fsolve(eval({e||(1..8)}, [y1[n]=iny1, y2[n]=iny2]), {vars}), tolerance = tolerance):

would have produced an error; for instance:

restart
dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0}:
abserr := 1e-6;
                            0.000001
dsol := dsolve(dsys, numeric, abserr=abserr):
Error, (in dsolve/numeric/an_args/SC) keyword was '0.1e-5', optional keyword must be one of 'abserr', 
'delaymax', 'delaypts', 'differential', 'event_doublecross', 'event_initial', 'event_iterate', 
'event_maxiter', 'event_pre', 'event_project', 'event_relrange', 'event_stepreduction', 'events', 
'implicit', 'initstep', 'interpolate', 'interr', 'maxfun', 'maxstep', 'minstep', 'optimize', 
'output', 'projection', 'range', 'relerr', 'startinit', 'steppast'

To avaoid this error the correct syntax is (note the simple quotes arround the first occurence of abserr)

dsol := dsolve(dsys, numeric, 'abserr'=abserr):

Both the rhs have potentially 2 singularities: one at R=0, the other at R=1.
Computing limit(rhs(Sys[2]), R=1) returns alpha^2 and thus there is no problem at R=1 for the Phi-ode.

To resolve the 3 remaining singularities I propose to replace localy the rhs by its series expansion.
More precisely:

• T-ode:
  • replace the rhs at R=0 by a series expansion T_SL in the interval [0, h],
  • note T_CS the rhs in the interval [h, 1-h],
  • replace the rhas at R=0 by a series expansion T_SR in the interval [1-h, L), where L is any number > 1,
    (possibly L = 1+h and a fourth range [1+h, L) could be included where T_CS2 is equal to the rhs)
  • dsolve exactly {diff(T(R), R) = T_LS, T(0)=0.5},
  • dsolve numerically {diff(T(R), R) = T_CS, T(h)=U} where U is the value of T at point R=h for the previous solution,
  • dsolve exactly {diff(T(R), R) = T_RS, T(1-h)=V} where U is the value of T at point R=h for the previous solution.
    (possibly dsolve numerically in the fourth range)
• Phi-ode
  • replace the rhs at R=0 by a series expansion Phi_SL in the interval [0, h],
  • note T_CS the rhs in the interval [h, L), where L is any number > 1,
  • solve exactly {diff(T(R), R) = T_LS, T(0)=0.5},
  • solve numerically {diff(T(R), R) = T_CS, T(h)=U} where U is the value of T at point R=h for the previous solution.

Merge the pieces of solution.
I took arbitrary h=0.25, but this value can be adjusted; a comparison with @Preben Alsholm's solution is also presented
Finally, I've chosen alpha=1, T(0)=1/2 and Phi(0)=0.
 

dsolve_series_1.mw

Assuming both alpha and k are strictly positive:

Jn := (alpha, k) -> k*Int((BesselK(0, k*Zeta)/Zeta^3)*Zeta*BesselK(1, alpha*Zeta), Zeta = 1 .. infinity):
r  := rand(0. .. 10):
ak := r(), r();
evalf(Jn(ak))
                    9.631107879, 5.713953505
                                      -8
                        1.463786783 10  

Multiply this value by 4*h1*lambda to get the final result

You can define cot(x) as 1/tan(x) or cot(x) as tan(Pi/2-x) and even cot(x) as -tan(x+Pi/2) (and while not say that cot(x) = cos(x)/sin(x) ?).
In your case (where does the value 90 come from?) you can obtain evry different simplifications: they all depend on the rules you apply and the order you apply them.

Result y is, IMO, even simpler and easier to read than MMA's result. It uses the subjective simplification tan = 1/cot.

Result z is closer to MMA's but replaces Pi/2 by Phi in order to avoid earlier (and autimatic) simplifications. It uses the subjective simplification cot = 1/tan (which is just one of the many definitions of cot).

The last result starts from z and uses nother subjective simplification cot(Pi/2+alpha/2)=-cot(Pi/2-alpha/2).

restart:
x := -sin(alpha)*(L + (e + r)/tan(-Pi/2 + alpha/2))/2 + sin(alpha)*L:
eval(%, cot=(u -> 1/tan(u))):
y := collect(%, sin);
           /1     1            /1      \\           
           |- L + - (e + r) tan|- alpha|| sin(alpha)
           \2     2            \2      //           

x := -sin(alpha)*(L + (e + r)/tan(Phi + alpha/2))/2 + sin(alpha)*L:

eval(%, tan=(u -> 1/cot(u))):
z := collect(%, sin);
        /1     1            /      1      \\           
        |- L - - (e + r) cot|Phi + - alpha|| sin(alpha)
        \2     2            \      2      //           
eval(z, cot(Phi+alpha/2)=-cot(Phi-alpha/2));
        /1     1            /      1      \\           
        |- L + - (e + r) cot|Phi - - alpha|| sin(alpha)
        \2     2            \      2      //           

simplify.mw

s := singular(1/sin(x), x):
var := (indets(%) minus {x})[]:
eval(s, var = `#mo("ℤ")`);

As @tomleslie said, one BC is missing.

An example with a periodic BC for D[1](u) ( = diff(u(x, t), x) )

restart:
PDE  := diff(u(x, t), t) = -u(x, t)*diff(u(x, t), x) + 0.1 * diff(u(x, t), x$2):
IC   := u(x, 0) = sin(x):
BC   := u(0, t)=u(2*Pi, t), D[1](u)(0, t)=D[1](u)(2*Pi, t):
IBC  := {IC, BC}:

pds := pdsolve(PDE, IBC, numeric, time=t, range=0..2*Pi):  # just do as explained in the help page
plots:-display(
  seq(
    pds:-plot(
      t=tau
      , numpoints=50
      , color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
      , legend=('t'=nprintf("%1.2e", tau))
    )
    , tau in [seq](0..2, 0.25)
  )
)

Explanation

P := plot3d( [ b, ( b^2 ) / 4, -b / 2 ], b = -4 .. 4, c = -4 .. 4, thickness = 5, color = red  ); 

P is made of two elements

op(P) = MESH(...., COLOR(....)), THICKNESS(...)

When a MESH structure is plotted the "construction lines" of the mesh are too (default choice).
To avoid plotting them (see PLOT3D help page) you have to use the option STYLE(HIDDEN).

Then (no need to the THICKNESS option)

PLOT3D( MESH( op(op(P)), STYLE(HIDDEN) ) )

does the job.

Not that simple but this explains why the curve you got is black: the black boundaries of each element of the MESH "overwrite" the red on their interiors.

First thing to do is to load (ExcelTools:-Import ; look the corresponding help page) your data 

Data := ExcelTools:-Import(....):

Determine the number of lines:

N := numelems(M[.., 1])

To split Data into a train base and a test base :

• define the fraction tau of data in the 
• set NT := round(tau*N)  (or floor or ceil instead)
• realise a random splitting:

  combinat:-randperm([$1..N]):
  TrainSet := %[1..NT]:
  TestSet  := %%[NT+1..N]

• then:

  TrainData := Data[TrainSet];
  TestData  := Data[TestSet];

Apply Statistics with the fit method (linear or nonlinear) to TrainData..
Evaluate the quality of the fit on TestData.

Provide an Excel data file if you want more details.

The reason why you get a wong picture for your second test case is detailed in the mw file I'm nor able to load right now (it will be the case in a few hours).


Concerning this second test case (Maple 2020)

restart:
with(plots):
with(ComputationalGeometry):

# Feasibility domain

Q := inequal(
       Or(u^2+4*v^2 <= 4, And(u^2+v^2 < 4, (u+2)^2+2*(v-1)^2 <= 2))
       , u=-2..2, v=-2..2
       , nolines
     ):

# Extract the polygons that <inequal> has built

POL := select(has, [op(Q)], POLYGONS):
numelems(POL): 
                         2

# Mapping & plotting

F := plottools:-transform((u, v)  -> [u^3-v^2, u^2-v^3]):
display(F(Q)):

# Searching noon at 2pm
# The analysis of the poorness of the above plot shows it comes from the closure of the drawing
# of some of the elemmentary polygons POL contains (this will be clear in the mw file).
# All these concerned polygons having more than 3 sides, the basic idea is to triangulate 
# the feasibility domain using DalaunayTriangulation.

# first component of POL

pol  := op(1..-4, op(1, POL)):
pts  :=  `<,>`(pol)
tris := DelaunayTriangulation(pts):
display(map(x -> plottools:-polygon(pts[x], style='polygon'), tris), axes=none):
d1 := display(F(%)):

# second component of POL

pol  := op(1..-4, op(2, POL)):
pts  :=  `<,>`(pol)
tris := DelaunayTriangulation(pts):
display(map(x -> plottools:-polygon(pts[x], style='polygon'), tris), axes=none):
d2 := display(F(%)):

# drawings

display(`<,>`(d1, d2, display(d1, d2))):

The result is already very good but could probably be improved if the construction of Q was refined.
Use optionsimplicit=[grid=[50, 50]] in inequal(...) ... but beware of the computation time for finer grids !!!


File to come (either from this logname or from my non professional logname "mmcdara")
Here it is

DrawingIssue.mw
Maple

Mathematica

Unless I'm mistaken (it's early in the morning and I haven't get my coffee) the answer is 260

Reasoning.

1. Choose a color for the right triangle : 5 choices.
2. Choose a color for top and bottom triangles : 4 choices for each.
3. For the left triangle :
   1. if colors of top and bottom triangles are the same, 4 choices are left.
   2. If they are not the same 3 choices are left.

The probability that the top and bottom triangles share the same color is 1/4.
Thus 
N = 5 * 4 * 4 * (1/4 * 4 + 3/4 * 3) = 5 * 4 * (4 + 9) = 20 * 13 = 260

Given the 2 elementary events R = "a Red ball is drawn" and W = "a White ball is drawn", this will give you all the compound events of an experiment in which 3 balls are drawn from the box.

restart:
alias(C=combinat):
Box := [R$3, W$3]:
(op @ C:-permute)~(C:-choose(Box, 3));

BTW, you code can be simply written this way

combinat:-powerset({$1..3})

As isolve seems (IMO) more powerful and versatile than msolve I have build a simple procedure _msolve  which uses syntax of msolve and the power of the isolve .
 

restart:
_msolve := proc(eqs, m)
  local v := convert(indets(eqs, name), list);
  local e := lhs~(eqs) =~ rhs~(eqs) +~ m *~ parse~(cat~(_K, [$1..numelems(v)]));
  local s := isolve(convert(e, set));
  map(u -> v =~ eval(v, u), [s]);
end proc:

_msolve(x^2=-1, 5);  # "your" problem
  [ [ x = 3 - 5 _Z1 ], [ x = 2 - 5 _Z1 ] ]

_msolve(2^i=3, 19);  # from msolve help page
  [ [ i = 13 - 18 _Z1 ] ]

_msolve({3*x-4*y=1, 7*x+y=2}, 19); # from msolve help page
  [ [ x = 15 + 57 _Z1 + 76 _Z2,  y = 11 + 38 _Z1 + 57 _Z2 ] ]

# A generalization of the previous system where 2 modulii are used:
  eqs := {3*x-4*y=1, 7*x+y=2}:
  m   := [7, 19]
_msolve(eqs, m);
  [ [ x = 15 + 37 _Z1 + 76 _Z2, y = 11 + 26 _Z1 + 57 _Z2 ] ]

# a specific instance
eval(%[], [_Z1=1, _Z2=2]);
   [ x = 204, y = 151 ]

# check
eval(irem~(lhs~(eqs), m), %);
   [ 1, 2 ]

Let P the list of (2D) points:

ch  := simplex:-convexhull(P, output=true):
pch := plottools:-getdata(ch)[1][1];    
n:=2:
while pch[n, 1] > pch[n_1, 1] do
  n := n+1:
end do:
NewtonPolygon := pch[1..n-1]

Explanation
pch contains the set of points (a matrix) ordered with respect to an increasing curvilinear abscissa, starting from the leftmost point and in counter clockwise sense.
The Newton polygon is the submatrix of pch wich contains all the points with increasing "x" (1st column of pch).
The loop is aimed to determine the nth point whose abscissa "x" is lower than the previous one

A minimalist answer you can decorate as you want.
I hope there is no typos for I do not have Maple right now.

(If you have any problems I will reply you when I'm home [login = mmcdara])

restart:
with(plots):
with(Maplets[Elements]):

f := proc(p1, p2)
local maplet:
maplet := Maplet(
  [
    [
      Plotter['PL1'](p1),
      Plotter['PL2'](p2)
    ]
  ]
):

Maplets[Display](maplet):
end proc:

A := dualaxisplot(sin(x), cos(x), x=0..2*Pi):
B := dualaxisplot(sin(x), cos(x), x=0..2*Pi, color=[green, blue]):

f(A, B);