@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

@Alejandro Jakubi  Your wonderful programmatic way  works, for example, in Standard Maple 12 on Windows 8.1:

cat(`&#`, convert("25AA", decimal, hex), `;`);

@Markiyan Hirnyk   Without a doubt, your plot out of competition!

@Rouben Rostamian   The parentheses are necessary.

Test the solution:

P:=Matrix([[ 0 , .5 , .5 , 0 , 0 , 0 ], [ 1/3 , 0 , 0 , 1/3 , 1/3 , 0 ], [ 1/3 , 0 , 0 , 0 , 1/3 , 1/3 ], [ 0 , 1 , 0 , 0 , 0 , 0 ], [ 0 , .5 , .5 , 0 , 0 , 0 ], [ 0 , 0 , 1 , 0 , 0 , 0 ]]):

pii:=Vector[row]([ a , b , c , d , e , f ]):

solve({seq((pii.P)[i]=pii[i], i=1..6)});

eval([seq((pii.P)[i]=pii[i], i=1..6)], %);

                                                       {a = 2.*f, b = 3.*f, c = 3.*f, d = f, e = 2.*f, f = f}

             [2.000000000*f = 2.*f, 3.0*f = 3.*f, 3.0*f = 3.*f, 1.000000000*f = f, 2.000000000*f = 2.*f, 1.000000000*f = f]