otat := proc(a::`=`, b::list)

     local q;

     q := isolate(a, y);

seq(eval(q, b[i]), i=1..nops(b));   # returns values of  y  for all values of  x  in the list  b    

end proc:

Example:

otat(x-y = 1, [x = 10, x = 11, x=14]);

                  y = 9, y = 10, y = 13

Addition:  The line  seq(eval(q, b[i]), i=1..nops(b));   can be written shorter  eval~(q, b)[];   with the same output

Boundary problems can rarely be solved symbolically. And your task even has  a parameter K. Such problems are solved numerically for a particular value of the parameter:

K:=1:

Eq := diff(y(x), x, x) = -(x^2+1)*y(x)+K;

sol:=dsolve({Eq, y(-1) = 0, y(1) = 0}, y(x),numeric);

plots[odeplot](sol,[x,y(x)], x=-1..1, color=red, thickness=2);

Here are another 6 ways (from the shortest to the longest) for the same list  L:=[7, 9, 13, evalf(Pi), a, b]:

L^~2;

`^`~(L,2);

map(`^`,L,2);

[seq(t^2, t=L)];

[L[i]^2 $ i=1..nops(L)];

[seq(L[i]^2, i=1..nops(L))];

restart:

with(plots):

f:=2*x:

g:=x^2: 

h:=unapply(piecewise( 0<=x and x<= 5, f, 5<x and x <= 10, g ), x):

A:=[1,2,8]:

h_table_A:=map(h,A);

                                         h_table_A := [2, 4, 64]

Addition: another way (without  unapply)

restart:

with(plots):

f:=x->2*x:

g:=x->x^2:

h:=x->piecewise( 0<=x and x<= 5, f(x), 5<x and x <= 10, g(x) ) :

A:=[1,2,8]:

h_table_A:=map(h,A);

Or:

ans:=solve({alpha = (1/2)*lambda*alpha*Pi, beta = (1/2)*lambda*beta*Pi}, lambda);

assign(ans);

lambda^6,  lambda^8;

local Prime:

Prime:= proc (a, b)

local i, p, r;

p := {};

r :={} ;

for i from a to b do

   if isprime(i) then p := p union {i} ;

     if type(sqrt(p[-1]-1), integer) then r := r union {p[-1]} end if; end if;

end do;

p, r;

end proc:

Your example:

Prime(1, 10);

                       {2, 3, 5, 7}, {2, 5}

Another example:

Prime(1, 100);

   {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}, {2, 5, 17, 37}

I guess that you are interested only real solution:

RealDomain[solve](convert({0.001=x*exp(0.6/(n*0.026)), 0.01=x*exp(0.68/(n*0.026))}, fraction),{x,n});

evalf(%);

n[1] := 1.5:  n[2] := 1:

sin(theta[g]):='n[2]'/'n[1]' = cat(convert(n[2], symbol)/convert(n[1], symbol), `=`, `0`, n[2]/n[1]);

Also, the initial value for  P  should be  0  instead of  1.

I think it is more natural to use a formula instead of loops, for example as in Wiki

restart;

f:=x->2^x:

M:={-2,-1,0,1,2}: n:=nops(M):

L:=(i,x)->normal(mul((x-j)/(M[i]-j), j=M minus {M[i]})):

P:=add(f(M[i])*L(i,x), i=1..n);

plot([f(x), P], x=-2..2, color=[red,blue]);

We see that the approximation is very good (the graphics are almost identical).

You can prevent the calculation of  ln(1.5)  just treating  1.5  as a symbol:

diff(`1.5`^x,x);

y:=sin(x):  y1:=x*sin(x):

maximize(y, x=0..Pi);

evalf(maximize(y1, x=0..Pi));

Optimization[Maximize](y1, x=0..Pi);

                                             1

                                   1.819705741

        [1.81970574115965, [x = 2.02875783126017]]

eq1:=x^2-y^2=a*z^2:

eq2:=x-y=a*z:

factor(eq1/eq2);

              x + y = z

S:={seq(seq(i*j, i=7..9997,10), j=1..9991,10)}:  # The set of all the products in increasing order

S[1..10];  # The first 10 products

S[-10..-1];  # The last 10 products

                                                 {7, 17, 27, 37, 47, 57, 67, 77, 87, 97}

{99580117, 99580297, 99580377, 99580437, 99680087, 99680207, 99680247, 99780057, 99780117, 99880027}

B:=[A1,A10,A11,A12,A2,A3,A4,A5,A6,A7,A8,A9]:

sort(B,(x,y)->parse(convert(x,string)[2..-1])<parse(convert(y,string)[2..-1]));

                 [A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12]

restart;

w:=1+t/100:

sol:=dsolve({diff(x(t),t$2)+w^2*x(t)=0, y(t)=diff(x(t),t), E(t)=0.5*(y(t)^2+w^2*x(t)^2), J(t)=E(t)/w, x(0)=1, D(x)(0)=2}, {x(t),y(t),E(t),J(t)}, numeric):

 plots[odeplot](sol, [[t,E(t)], [t,J(t)]], t=0..10, color=[red,blue], thickness=2, legend=['E(t)','I(t)']);