@David Sycamore 

What could possibly be interesting about this process and the code that performs it--specifically the part about checking for solutions up to a finite limit--given that I have proven beyond any doubt that for every case where there has been apparently no solution that indeed there is no solution---no matter how high the limit would be raised? Do you not understand my proof, or believe it to be incorrect?

This comment is directed at David, not Tom, who was merely providing the code to David's specifications. It's those specifications that now seem trivial or facile in light of what's been discussed in this thread.

@SirFrancisBacon 

Please use exactly the commands shown in ThU's Answer above. 

V:= D__e*(1 - exp(-a*(R - R__e)))^2;  #not V(R):= ....
limit(V, R= infinity) assuming a>0;    #not V(R)

@AndrewG LPSolve will only return one solution, even if there are more than one; however, I wouldn't say that it does so at random.

If you did want to list the feasible set, then seq (as shown by Kitonum) is a good way. But note that there are many such problems where that would be impossible due to the size of the set, yet for which LPSolve will provide the solution.

@Christian Wolinski The way that igcdex does it is cheating IMO. It first gets the GCD by calling builtin igcd, and then it backsolves for the Bezout coefficients by computing a modular inverse^[*1]. In the sense of analysis of algorithms, it's more efficient to accumulate the Bezout coefficients during each step of the Euclidean algorithm. However, the way that igcdex does it does have a substantially better average-case time in Maple simply due to igcd being builtin and part of the superfast GMP package. 

If GMP doesn't have an igcdex, that's a bit surprising; and if it does, then why isn't it used in Maple?

[*1] ... whereas the standard algorithm for a modular inverse uses the Bezout coefficients. 

@David Sycamore There may be independent copies which haven't been found yet, but the formula above only gives dependent copies for n > 1. In other words, consider these 3 sequences:

1. a(q) $ q= 0..infinity  (n=0, the primary sequence)
2. a(4*(q+2)-2) $ q= 0..infinity  (n=1, 1st-level identical subsequence)
3. a(16*(q+2)-2) $ q= 0..infinity  (n=2, 2nd-level identical subsequence).

Now #3 is a subsequence of #2 (as well as being a subsequence of #1) with exactly the index pattern relative to #2 as #2 has as a subsequence of #1. (These formulas are all valid for q=0, btw.)

Since my algebraic formulation is not mathematically proven to be equivalent to the sequence's definition (and I agree with VV that a proof would be difficult), you should post both versions (table-based and algebraic) of my Maple code on OEIS. Here are my polished copies that you may use:

a:= module()
description
     "Table-based formulation of David Sycamore's selfie sequence---"
     "straight from its definition"
;
option
     `Author: Carl J Love, 13-Oct-2019`
;
export
     #records index (position) where values occur, other than most recent
     pos:= table([0= 0]),
     nextpos:= table([0= 1]), #most-recent position for each value
     used:= table(sparse, [0= 2]), #times value has been used
     ModuleApply:= proc(n::nonnegint)
     option remember;
     local
          #`thisproc` is a keyword for recursive calls. 
          p:= thisproc(n-1), 
          r:= `if`(used[p] > 1, pos[p], thisproc(p-1))
     ;
          #shuffle down:
          (pos[r], used[r], nextpos[r]):= (nextpos[r], used[r]+1, n);
          return r
     end proc
;
     (ModuleApply(0), ModuleApply(1)):= (0,0) #Initialize remember table.
end module
:#------------------------------------------------------------------------
A:= proc(n::nonnegint)
description
     "Fast algebraic formulation of David Sycamore's selfie sequence, "
     "based on empirical arithmetic analysis---not mathematical proof. "
     "Works for very large n."
;
option 
     `Author: Carl J Love, 13-Oct-2019`, 
     remember
;
local 
     N:= 2^(ilog2(n+1)-2), #greatest power of 2 not exceeding (n+1)/4
     r, q:= iquo(n,4,r)-1, #q = (n div 4) - 1, r = n mod 4
     p, h:= iquo(r,2,p)+1  #p = if(n::even, 0, 1), h = if(r=1, 1, 2) 
;
     #For integer x, sign(x) = 2*Heaviside(x)-1 = if(x < 0, -1, 1).
     #The `thisproc` is a keyword for a recursive call.
     `if`(p=0, `if`(r=0, q, thisproc(q)), (6-sign(6*N-n-h))*N-2*q-4-h)
end proc: 
#Initialize remember table:
(A(0), A(1), A(2), A(5)):= (0, 0, 0, 2)
: 

@ANANDMUNAGALA To see the objective and constraints in algebraic form, simply enter this at any command prompt after executing the worksheet that I posted:

<P.V = max, C.V <=~ R>;

@Kitonum Thank you, that's good to know. It looks like the simplex package allows similar convenient syntax.

@David Sycamore Any sequence which contains itself as a proper subsequence necessarily contains infinitely many copies of itself. Thus, IMO, this "contains infinitely many copies of itself" aspect is much too trivial to mention in any formal article or OEIS entry.

Regarding the iterated or "generalised" spacing: It goes by powers of 4, not 2. We have

a(q) = a(4^n*(q+2) - 2), q >= 0, n >= 0.

@David Sycamore Thanks for the clarification, especially of last. My first several attempts at generating the sequence failed because I interpretted "last" as "most recent" rather than "most recent unmarked."

@David Sycamore

Your new sequence (the one generated by the code in my Answer immediately above) contains a copy of itself as an evenly spaced subsequence; specifically, a(4*q+2) = a(q-1), q > 0. I derived simple algebraic formulae for the other remainders mod 4: a(4*q+r), r in {0, 1, 3}. Thus, the sequence can be computed by an ordinary recursion (no counting, auxilliary tables, etc.), and recursion is only used in 1/4 of the cases.

A_:= proc(n::nonnegint)
option remember;
local N:= 2^(ilog2(n+1)-1), r, q:= iquo(n,4,r)-1;
   `if`(r=0, q, `if`(r=2, thisproc(q), `if`(n > 3*N-r/2, 7, 5)*N/2 - 2*q - (r+9)/2))
end proc: 
#Initialize remember table:
(A_(0), A_(1), A_(2), A_(5)):= (0, 0, 0, 2)
: 

I verified that this matches the previous code for the first 100,000 terms. This code computes about 300,000 terms per second on my modest computer.

This sequence needs a catchy name.

@ajfriedlan Issue the command currentdir(); before both the save and read commands. Do the two currentdirs return the same string? If not, that's your problem. Use currentdir("some path name") to set the current directory to the same value in each worksheet.

@tomleslie I'm not sure---and I wouldn't mind being proven wrong---but I don't think that that'll work with a Matrix of highly symbolic entries such as described by the OP.

@ajfriedlan Sorry, I missed a comma in my command. It should be

save M, "MyData.map";

@David Sycamore I need you to clarify the rule "a(n+1) = index of the last term k = a(n)". In particular, what is meant by index and what is meant by last? I can think of two conflicting definitions for each:

Index: For the term a(n), is n the index? Or is n+1 the index because the term appears in the (n+1)^th ordinal position?

Last: In the sequence 1, 2, 3, is 1 or 3 the last term?

@Kitonum I hope that you don't need me to tell you that that's extremely inefficient; for example, it takes over 500 times as long as mine to compute the first 2^13 terms.