@Carl Love Thank you, the problem is solved.

@Carl Love I have 

n := 3:
i := 1:
j := 2:
T[-3] := 0: T[-2] := 0: T[-1] := 0: T[0] := 0:
T[1] := 0.1: T[2] := 0.25: T[3] := 0.35: T[4] := 0.65:
T[5] := 1: T[6] := 1: T[7] := 1: T[8] := 1:

print(int(t-> N(i, n, t, T)*N(j, n, t, T), 0..1, numeric));

Maple returns int(t -> N(i, n, t, T) N(j, n, t, T), 0. .. 1.)

of course j is known as well

@Carl Love 

The reverse order didn't solve the problem. 

p1[i] are curves whereas pts[1], pts[2] are disks, i.e.,

pts1 := disk([P[2][1][1], P[2][1][2]], 0.009, color = pink): 
pts2 := disk([P[3][1][1], P[3][1][2]], 0.009, color = yellow):

p1[1] := plot( 
[
    add(B(n[1], i, x)*P[1][i+1][1], i = 0..n[1]),
    add(B(n[1], i, x)*P[1][i+1][2], i = 0..n[1]), 
    x = 0..1
],
    color=blue, linestyle=solid, axes=normal, scaling=constrained, numpoints=100, thickness=2
):

@Carl Love

I'm getting these equations while I'm trying to minimize a function of several variables. So probably, I'm going to use quadratic programming and nonlinear programming to find the minimum.

@Carl Love Thanks, I checked the fsolve procedure and, unfortunately, the results are not good enough.

Besides, I have also another system of equations (six equations of the third degree with six unknowns) and fsolve doesn't produce any solution at all.

@Markiyan Hirnyk To run DirectSearch-Search for the first time I need to pass the parameter initialpoint. How can I get it ?

@Markiyan Hirnyk To run DirectSearch-Search for the first time I need to pass the parameter initialpoint. How can I get it ?

Thanks. I've got one question. How do you get your first initialpoint ? Is it something that is random ?

Thanks. I've got one question. How do you get your first initialpoint ? Is it something that is random ?

@PatrickT 

EX := sum(
                      (
                          sum(
                                    PX[k+1]*B(n, k, T[j+1], c[1], c[2]),
                                    k = 0..n
                                )
                          -sum(
                                     WX[k+1]*B(m, k, T[j+1], c[1], c[2]),
                                     k = 0..m
                                 )
                      )^2,
                      j=0..nops(T)-1
                  ):

    SX := solve(
                   [
                       seq(diff(EX,WX[k+1])=0, k = 0..m), # these are linear
                       seq(AX[k+1]*(WX[k+1]-rightX)=0, k = 0..m), # these are nonlinear, but very simple
                       seq(BX[k+1]*(leftX-WX[k+1])=0, k = 0..m), # these are nonlinear, but very simple
                       seq(WX[k+1] <= rightX, k = 0..m),
                       seq(WX[k+1] >= leftX, k = 0..m),
                       seq(AX[k+1] >= 0, k = 0..m),
                       seq(BX[k+1] >= 0, k = 0..m)
                   ],
                   [
                       seq(WX[k+1], k = 0..m),
                       seq(AX[k+1], k = 0..m),
                       seq(BX[k+1], k = 0..m)
                   ]
               ):

where

n := nops(PX)-1:

PX := [221, 199, 123, 15, 97, 127, 219, 213, 186, 161, 213, 218]

c:= [0, 1]:

T := [seq(evalf(k/15), k = 0..15)]:

m := 7:

#Bernstein polynomials:

B := (n, i, x, t0, t1) -> binomial(n,i)*((x-t0)/(t1-t0))^i*(1-(x-t0)/(t1-t0))^(n-i):

leftX := 1:
rightX := 300:

WX, AX, BX are unknowns

So what I'm trying to do is to reduce the degree of Bezier curve with some constraints, minimize the least square distance between curves.

thanks, it works !

thanks, it works !