X:=[10,15,20,25,30,35,40]: 

Y:=[2.01,4.87,6.50,9.875,14.00,18.90,24.40]:

 f:=x->`if`(is(x in X), Y[ListTools[Search](x,X)], undefined):

Examples:

f(35);    f(1);

                  18.90

                undefined

If we assume that a solution exists and substitute into the equation  x = 0, we obtain  theta (0) = 0. But then, solving with these initial conditions, we have only the trivial solution:

eq:=(1+theta(x))*diff(theta(x),x)+x*diff(theta(x),x)^2+x*(1+theta(x))*diff(theta(x),x,x)-theta(x)^3-diff(theta(x),x)=0;

ics:=theta(0)=0, D(theta)(0) = 0;

Sol:=dsolve({eq, ics}, numeric);

plots[odeplot](Sol,[x,theta(x)], x=0..1, axes=box);

Here is step by step solution. I used my procedure  JordanGausse . This procedure step by step solves any system of linear equations by Gauss - Jordan elimination. In the procedure, the Russian text is used. After I prepare the English version, I will submit  code of the procedure to the mapleprimes forum.

Unfortunately, when loading the Russian text was distorted. Here is a translation into English:

Step 1: write the system in matrix form

Step 2: divide row 1 by 4

Step 3: row 1, multiplied by -2, add to row 2

Step 4: row 1, multiplied by -4, add to row 3

Step 5: swap the positions of columns 2 and 3

Step 6: divide row 2 by 3/2

and so on ...

Download Gauss-Jordan.mws

It is well known that a square system has a unique solution if and only if it's main determinant is not 0.

restart;

A:= < k+3, 2, k-4, 3;

      0,   2,  -9, 5;

      0,   0, k^2+k-2, k-1 >;

B:=A[..,1..3];

V:=A[..,4];

LinearAlgebra[Determinant](B);

Sol:=[solve(%)];

We have a unique solution for any number  k  not included in this list.

Next solve separately for each number from  list  Sol:

for i from 1 to 3 do

'k'=Sol[i],  LinearAlgebra[LinearSolve](op(eval([B,V],'k'=Sol[i])), free=C);

od;

Thus for  k=-3  and  k=1  we have infinitely many solutions, for  k=-2  no solutions

You made ​​a few syntax errors, which I corrected.
Your nonlinear system has 5 parameters. Maple can not solve it symbolically for all values ​​of the parameters, but Maple can  solve this system for the specified parameters.

restart;

f := proc (E2, T, DHT, Ca, Cshbg)

local KsT, KaT, KsE2, KaE2, KsDHT, KaDHT, eq1, eq2, eq3;

KsT := 0.10e11; KaT := 4.6*0.10e6; KsE2 := 3.14*0.10e10; KaE2 := 4.21*0.10e6; KsDHT := 3*0.10e6; KaDHT := 3.5*0.10e6;

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT));

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT));

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT));

solve({eq1, eq2, eq3})

end proc:

Example:

f(1, 2, 3, 4, 5);

We got 10 solutions.

Procedure  UniqueColumns  returns  the unique columns of a Matrix.

UniqueColumns:=proc(A::Matrix)

local L;

uses LinearAlgebra, ListTools;

L:=[Categorize((x,y)->Equal(x,y), [seq(A[..,i],i=1..ColumnDimension(A))])];

op(map(op,select(x->is(nops(x)=1), L)))

end;

Your example:

Q := Matrix([[1,1,1,0,0,0,-2,-1], [0,0,0,1,1,1,0,-3]]);

UniqueColumns(Q);

Write your function as a procedure:

f:=proc()

your expression

end;

Example:

f:=proc()

add(args[i]^2, i=1..nargs)

end;

f(1, 2, 3);

       14

I do not understand the meaning of the problem. If you want to write it in the same loop then 

for i in  [seq(i, i=0..100, 0.2), seq(i, i=100..0, -0.2)]  do ... od;

Multiplication signs, division sign and missing parenthesis

((5*x^2*y^2)-6*x^4)/((1-y^2)*(x^2+y^2)-(y^2*(1-x^2-y^2)));

The parenthesis can be put in another place.

This is arithmetical progression.

f:=unapply(rsolve({u(n+1)=u(n)+d, u(1)=a}, u(n))/b, n);

                              f := n -> (a+d*(n+1)-2*d)/b

Example:

eval(f(100000000000), {a=25, d=25, b=5});

                                         500000000000

Look at my posts   http://www.mapleprimes.com/users/Kitonum/posts 

You can find a number of interesting applications of them.

You can use geometric transformations from  plottools  package.

Examples:

with(plottools): with(plots):

F:=curve([[-4,0],[0,3],[1/2,0], [-4,0]], thickness=3):

F1:=reflect(F, [[-1,4],[1,4]]):

F2:=translate(F, -4, 2):

F3:=rotate(F, -Pi/2, [3,0]):

T:=textplot([[-1,1,"F"], [-1,7,"F1"], [-5,3,"F2"], [4,4,"F3"]], font=[TIMES,ROMAN,20], color=blue):

display(F, F1, F2, F3, T, scaling=constrained);

Code of procedure:

ExponentOfPower:=proc(monom, var)

`if`(not has(monom,var), 0,  `if`(nops(select(has,monom,var))>1, op(2,select(has,combine(monom),var)), 1));

end:

Examples:

term1:=3*x^k*y^(k+2): term2:=3*y^(k+2): term3:=3*x*y^(k+2): term4:=3*x^k*y^(k+2)*x^2:

op(map(ExponentOfPower, [term||(1..4)], x));

                                                                        k, 0, 1, k+2

Example with Carl's procedure:

power(x^x, x);

             x*ln(x)+x

f := t->piecewise(0 < t and t < 10, 1-t, 10 < t, t);

plot([f(t), [10, t, t = -9 .. 10]], t = 0 .. 15, thickness = [2, 1], color = [red, black], linestyle = [1, 2], discont, view = [-2.5 .. 14.5, -10 .. 14.5], scaling = constrained);

term:=3*x^k*y^(k+2):

op(2, select(has, term, x)),   op(2, select(has, term, y));

                                                  k, k+2