@Carl Love  I think that  discont  option  can only build one circle at the respective values  x. There is no need to build a hollow circle on the right end, because we consider the single-valued functions. For example, when  x = 2, the function value is 2, and the solid circle indicates it .

@Markiyan Hirnyk  Thanks. I think everyone can come up with many applications of the procedures  Composition  and  Composition1. Here is one of them:

The problem - how many ways 1 dollar (100 cents) can be exchanged  by 10 coins (in denominations of 1, 5, 10, 25, 50 cents) in any combination (exactly 10 coins should be used)?

Composition1(100,10, {1,5,10,25,50}):

convert(map(t->sort(t), %), set);

nops(%);

@Markiyan Hirnyk  for your brilliant solution!

@Markiyan Hirnyk  If my explanation is not enough for you, here's the direct check. For it, I added to the original equation 3 additional equations  diff(y(x), x$2)=k(x), diff(k(x), x)=p(x), diff(p(x), x)=q(x)

restart;

Digits:=20:

z:=x->diff(y(x),x,x)+sin(y(x)):

u:=diff(z(x),x):

v:=diff(z(x),x,x):

w:=diff(y(x),x$3)+diff(y(x),x$2):

sys:=subs({y(x)=a,diff(y(x),x)=b,diff(y(x),x$2)=c,diff(y(x),x$3)=d,diff(y(x),x$4)=0.1},{z(x)=0,u=0,v=0, w=-0.1});

fsolve(sys);

assign(%);

sol:=dsolve({z(x)=0,diff(y(x),x$2)=k(x),diff(k(x),x)=p(x),diff(p(x),x)=q(x),y(0)=a,D(y)(0)=b,k(0)=c,p(0)=d}, numeric);

sol(0);

%[2]+%[3], %[4];

@Markiyan Hirnyk  Let  y(x)  is a solution of the original equation, satisfying to conditions  diff(y(x), x$3)(0)+diff(y(x), x$2)(0)=-0.1  and  diff(y(x), x$4)(0)=0.1

Substitute  y(x) into the original equation and differentiate it 2 times. We got 3 identities that hold for any x, in particular, for x=0. Introduce the notations for  y(0)  ( a=y(0) ) and for all derivatives at  x=0  and adding your initial condition, we obtain a system of 4 equations with 4 unknowns. The further solution is obvious.

You wrote "Is it possible to verify those conditions by the fdiff command?" Directly is impossible, because  y(x)  should be an algebraic expression. A possible approach - using individual points  [x, y(x)]  firstly to get this expression, for example, by interpolation.

@Markiyan Hirnyk  These equations are your initial conditions (to calculate the arbitrary constants   _C1  and  _C2  in symbolic solution  sol  in my first post in this thread).  Also pay attention to the addition to my second post.

Your example has no solutions, because we have inconsistent system:

restart;

sol:=rhs(dsolve(diff(y(x),x,x)+sin(x)=0));

{eval(diff(sol,x$3),x=0)+eval(diff(sol,x$2), x=0)=-0.1, eval(diff(sol,x$4), x=0)=0.1};

                                               sol := sin(x)+_C1*x+_C2

                                                   {-1 = -.1, 0 = .1}

@toandhsp   Here are the solutions with  Composition1  procedure:

N1:=0:

for n from 4 to 296 by 2 do

N1:=N1+nops(select(t->nops(t)=nops(convert(t,set)), Composition1(n,3, {$ 1..100})))/3!;

od:

N1;   # Solution of the problem 1

N2:=0:

for n from 6 to 15 do

N2:=N2+nops(select(t->nops(t)=nops(convert(t,set)) and irem(`+`(op(t)),3)=0, Composition1(n,3, {$ 1..6})));

od:

N2;   # Solution of the problem 2

Of course, these examples are easier and faster to solve in another way:

binomial(50,3)+50*binomial(50,2);   # Solution of the problem 1

select(t->irem(`+`(op(t)),3)=0, combinat[permute](6,3)); 

nops(%);   # Solution of the problem 2

@Carl Love  I meant to set up an appropriate graphical structure (plot or plots[implicitplot]) with options like grid  or  numpoints  so  that independent variable were passing with predefined step, for example,  t=-0.05, -0.049, -0.048, ... .

@Carl Love  When we extract the data from a graphical structure, can we make so that they were passing with predefined step, for example,  t=-0.05, -0.049, -0.048, ... ?

@acer  Good idea!

Composition1  procedure improves the previous  Composition  procedure. Added option  S  - set, which includes elements of the composition. By default  S = {$ 0.. n} .  In this case  Composition1(n, k, res)  is equivalent to  Composition(n, k, res) .

Composition1 := proc (n::nonnegint, k::posint, res::{range, nonnegint} := 0, S::set:={$0..n})

local a, b, It, L0;

if res::nonnegint then a := res; b := n-(k-1)*a  else a := lhs(res); b := rhs(res) fi;

if b < a or b*k < n then return `No solutions` fi;

It := proc (L)

local m, j, P, R, i, N;

m := nops(L[1]); j := k-m; N := 0;

for i to nops(L) do

R := n-`+`(op(L[i]));

if R <= b*j and a*j <= R then N := N+1;

P[N] := [seq([op(L[i]), s], s = {$ max(a, R-b*(j-1)) .. min(R, b)} intersect S)] fi;

od;

[seq(op(P[s]), s = 1 .. N)];

end proc;

L0 := [[]];

(It@@k)(L0);

end proc:

Two examples:

Composition1(20,3, {seq(2*i, i=1..10)});  # Elements should be posint and even

[[2, 2, 16], [2, 4, 14], [2, 6, 12], [2, 8, 10], [2, 10, 8], [2, 12, 6], [2, 14, 4], [2, 16, 2], [4, 2, 14], [4, 4, 12], [4, 6, 10], [4, 8, 8], [4, 10, 6], [4, 12, 4], [4, 14, 2], [6, 2, 12], [6, 4, 10], [6, 6, 8], [6, 8, 6], [6, 10, 4], [6, 12, 2], [8, 2, 10], [8, 4, 8], [8, 6, 6], [8, 8, 4], [8, 10, 2], [10, 2, 8], [10, 4, 6], [10, 6, 4], [10, 8, 2], [12, 2, 6], [12, 4, 4], [12, 6, 2], [14, 2, 4], [14, 4, 2], [16, 2, 2]]

Composition1(30,4, select(isprime,{$3..30}));  # Elements should be prime greater 2

[[3, 3, 5, 19], [3, 3, 7, 17], [3, 3, 11, 13], [3, 3, 13, 11], [3, 3, 17, 7], [3, 3, 19, 5], [3, 5, 3, 19], [3, 5, 5, 17], [3, 5, 11, 11], [3, 5, 17, 5], [3, 5, 19, 3], [3, 7, 3, 17], [3, 7, 7, 13], [3, 7, 13, 7], [3, 7, 17, 3], [3, 11, 3, 13], [3, 11, 5, 11], [3, 11, 11, 5], [3, 11, 13, 3], [3, 13, 3, 11], [3, 13, 7, 7], [3, 13, 11, 3], [3, 17, 3, 7], [3, 17, 5, 5], [3, 17, 7, 3], [3, 19, 3, 5], [3, 19, 5, 3], [5, 3, 3, 19], [5, 3, 5, 17], [5, 3, 11, 11], [5, 3, 17, 5], [5, 3, 19, 3], [5, 5, 3, 17], [5, 5, 7, 13], [5, 5, 13, 7], [5, 5, 17, 3], [5, 7, 5, 13], [5, 7, 7, 11], [5, 7, 11, 7], [5, 7, 13, 5], [5, 11, 3, 11], [5, 11, 7, 7], [5, 11, 11, 3], [5, 13, 5, 7], [5, 13, 7, 5], [5, 17, 3, 5], [5, 17, 5, 3], [5, 19, 3, 3], [7, 3, 3, 17], [7, 3, 7, 13], [7, 3, 13, 7], [7, 3, 17, 3], [7, 5, 5, 13], [7, 5, 7, 11], [7, 5, 11, 7], [7, 5, 13, 5], [7, 7, 3, 13], [7, 7, 5, 11], [7, 7, 11, 5], [7, 7, 13, 3], [7, 11, 5, 7], [7, 11, 7, 5], [7, 13, 3, 7], [7, 13, 5, 5], [7, 13, 7, 3], [7, 17, 3, 3], [11, 3, 3, 13], [11, 3, 5, 11], [11, 3, 11, 5], [11, 3, 13, 3], [11, 5, 3, 11], [11, 5, 7, 7], [11, 5, 11, 3], [11, 7, 5, 7], [11, 7, 7, 5], [11, 11, 3, 5], [11, 11, 5, 3], [11, 13, 3, 3], [13, 3, 3, 11], [13, 3, 7, 7], [13, 3, 11, 3], [13, 5, 5, 7], [13, 5, 7, 5], [13, 7, 3, 7], [13, 7, 5, 5], [13, 7, 7, 3], [13, 11, 3, 3], [17, 3, 3, 7], [17, 3, 5, 5], [17, 3, 7, 3], [17, 5, 3, 5], [17, 5, 5, 3], [17, 7, 3, 3], [19, 3, 3, 5], [19, 3, 5, 3], [19, 5, 3, 3]]

@brian bovril  I did not include the number  0  as a valid integer in heading of  the procedure  NumbersGame  because with this number sometimes problems arise (though very rarely). If you replace the heading of  the procedure with

NumbersGame:=proc(Result::{integer,fraction}, Numbers::list(nonnegint), Operators::list:=["+","-","*","/"], NumbersOrder::string:="strict order", Parentheses::symbol:=no)

then all your examples easy to solve even without changing the order of the numbers  in the list  [0,1,2,5] . Only the representation of the number 1 requires a change order:

for n from 1 to 10 do

NumbersGame(n, [0,1,2,5]);

od; ``;

[NumbersGame(1, [0,1,2,5], "arbitrary order")][-7];  # Representation of 1

@Carl Love  Your code is compact and elegant, so vote up. But, unfortunately, it does not work in older versions.

In Maple 12:

P:= (n::posint, m::posint)->

     (x-> (d-> `if`(d[1]=d[-1], x, ``))(convert(x, base, 10)))~

          (LinearAlgebra:-RandomMatrix(n,m, generator= rand(100..999))):

          Error, missing operator or `;

@smith_alpha   1) If you just need to solve a system of equations, I do not think there are any advantages of PolynomialSystem  command comparing with  solve  command. Here is a quote from help on  SolveTools  package 

"The SolveTools package is a collection of building blocks for solving systems of algebraic equations. The routines from this package are used at the heart of the Maple solve command. An expert user can take advantage of the SolveTools package in order to perform individual steps towards solving a system of algebraic equations, thus allowing more control over how a solution is found".

2)  %[1]  is the first solution of the system, ie,  {x = 2, y = 1}

     eval([x, y], %[1])  is the list  [2, 1]

     eval([x, y], %[1])[]   just removes the square brackets, ie, we obtain a sequence  2, 1 .  This is necessary

     for     multiple assignment    a, b := 2, 1