The units of variance are always the squares of the units of the mean (expected return). So, using your current constraint budget = 10000, the expected return is on the order of a thousand dollars or so, and so the variance is on the order of a million or so. So, your constraint variance <= 0.02 is totally unrealistic and unachievable. You are thinking in terms of squared coeffient of variation (variance/mean^2) for the right side of that constraint whereas the left-side computation is the absolute variance. The former is unitless.

Two solutions:

Use unit budget (b:= 1): If you make the budget = 1, then the solution values for x will represent the fraction of your total budget to invest in each asset, and the solution for the objective function will be the growth rate.

Use coefficient of variation: Do this:

relvar:= 1.; #target relative variance
variance:= expand(X^%T.Q.X - relvar*growth^2) <= 0;

Note that this constraint is equivalent to sigma/abs(mu) <= sqrt(relvar) = relative standard deviation. I think that a target relative standard deviation of sqrt(0.02) = 14% shows a high aversion to risk, but that's a personal decision. Experiment with different values of relvar. In either case, DirectSearch:-Search seems to work better than Optimization:-NLPSolve.

This is a formula that all computer programmers should memorize: For integers a and b, the inclusive number of integers between and including them is abs(b-a) + 1. I call it the fence-post formula because it can be used to answer this: How many fence posts are needed to build a 10-meter fence with posts every 1 meter? Answer: 11. The failure to apply this formula correctly is one of the most common causes of bugs in computer programs. 

In most cases where there is apparently no solution, it is easy to prove that there is indeed no solution.

Proposition: Let p be a prime. Let e[0]=0, and let e[1], ..., e[p-1] be p-1 distinct positive integers such that  p + e[i] is prime for 0 <= i < p. Compute all the remainders e[i] mod p, and suppose that this set of remainders is {0, 1, ..., p-1}. Let q > p be any integer (prime or otherwise). Then the set Q = {q + e[i] | 0 <= i < p} necessarily contains a multiple of p. So, there's a composite in Q.  I urge you as a student of number theory to prove this proposition.

Now let's apply this proposition to the specific cases that you found, 105, 195, 231: The following computations are trivial to do by hand or even mental computation, but let's use Maple:

for n in [105, 195, 231] do
   P:= numtheory:-factorset(n);
   p:= min(P);
   E:= P -~ p;
   R:= irem~(E,p);
   if nops(R) = p then
      printf("Provably no solution for n = %d", n)
   fi
od;
                         P := {3, 5, 7}
                             p := 3
                         E := {0, 2, 4}
                         R := {0, 1, 2}
Provably no solution for n = 105
                        P := {3, 5, 13}
                             p := 3
                        E := {0, 2, 10}
                         R := {0, 1, 2}
Provably no solution for n = 195
                        P := {3, 7, 11}
                             p := 3
                         E := {0, 4, 8}
                         R := {0, 1, 2}
Provably no solution for n = 231

Now, it's probably obvious that for "small" n for which there is no solution, p is most likely to be 3. But it's not always 3. I believe that the first n for which it's not 3 is 95095, for which it's 5. Put 95095 into the for-loop above. Then P = {5, 7, 11, 13, 19}, and for every positive integer e, Q:= P +~ e necessarily contains a composite multiple of 5. Q may also contain a multiple of 3, but for e=6, it does not. 

Given a prime p, the following procedure finds a "no-solution" example whose smallest prime factor is p:

#Find a "no-solution" example with a given smallest prime, p:
FindExample:= proc(p::prime)
local R:= table([0= NULL]), n:= p, p1:= p, r:= 0, k;
   for k to p-1 do
      while assigned(R[r]) do 
         p1:= nextprime(p1);
         r:= irem(p1,p)
      od;
      R[r]:= NULL;
      n:= n*p1
   od;
   n
end proc
:
FindExample~([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]);
[ 6, 105, 95095, 215656441, 5070519693819151, 121621899527020215557, 
  21664759328930363899339259137, 6118163487924060025029213220329793, 
  3042203271558217610793849843552106188681707, 
  410050702193824620789441168165116752238670790985538614407, 
  243521388535217606004348248409375049907078159488916219637351331, 
  65167879497494548612895080562902666093018940639437608943583983945258079733989
  ]

For each of these, there's a good chance that the given example is the smallest such example, but I won't guarantee it.

The values in Q and r are empirical, by which I mean that they come directly from real-world observation (or are made up as a academic example) rather than calculation. Certainly, r is not calculated in the worksheet. The r-values are the expected return rates of the three assets (x, y, z) in the proposed portfolio.

The matrix Q is covariances, not returns (or return rates). It represents that the assets do not grow independently of one another. For example, the negative value (-.20) in the [1,2] position indicates that asset 1 tends to go up when asset 2 goes down and vice versa (covariance matrices are always symmetric).

The expected return (what you're calling EV) is X.r (vector dot product). The author calculates this as add(X[i].r[i], i= 1..3), which is unnecessary because that's exactly the same as X.r.

Your "average" method for EV is invalid (unless perchance the same amount of each asset is purchased, in other words x=y=z).

There is no way to compute r from Q---the concept is not meaningful. It would be analogous to computing the mean of a data set only knowing its standard deviation.

I think that you may have mistaken the letter Q as representing something measured quarterly (i.e., every three months). I don't know why that letter was used, but it has nothing to do with calendar quarters. 

All pictures are rectangles. Your picture is not just a penny; it's a penny within a black square background. That entire square is being squeezed into your plot of a circle. The way around this is to plot an entire square with the values outside the circle being undefined: 

plot3d(
   (x,y)-> `if`(x^2+y^2 <= 1, 0, undefined), -1..1, -1..1, 
   image= "penny.png", glossiness= 0, lightmodel= light4
)

This method is much simpler and much more efficient than those presented in the other two Answers (Kitonum and Tom Leslie):

i:= 1 + ListTools:-BinaryPlace(M[.., j], x, `<=`)

where M is the matrix, j is its sorted column, and x is the target value.

The command is

Sols:= solve(x^4-2*x^3-x^2+8*x+8-2*t, x, parametric, real)

This shows 0 solutions for t < 1, 1 for t=1, and 2 for t > 1. You can get Maple to express that directly by

subsindets[flat](Sols, list, nops);

The worksheet below has a module with a thorough implementation of my algorithm to find all (unique) expressions of a nonnegative integer n as a sum of four squares. It can also be used to get a specified number of solutions.

The key ideas of the algorithm:

1. Jacobi's four-square theorem makes it remarkably simple to count the total number of solutions (including multiplicities due to negatives and permutations).
2. The multiplicity of a 4-tuple is easy to count.
3. Generate solutions at random until the total is reached.

All details and references (all online) are in the comments in the code:

(*>>>----------------------------------------------------------------------
 * A module to compute expressions of an integer as a sum of four squares *
 * -----------------------------------------------------------------------*
 *														    *
 * This code requires Maple 2019 or later syntax.					    *
 *---------------------------------------------------------------------<<<*) 
Sums4Squares:= module()
description `Expressions of an integer as a sum of four squares`;
option
	`Author: Carl Love <carl.j.love@gmail.com> 24-Sep-2019`,
	`References: `,
		`"Lagrange's four-square theorem"`,								#Ref 1
			`https://en.wikipedia.org/wiki/Lagrange%27s_four-square_theorem`,
		`"Sum of two squares theorem"`, 									#Ref 2
			`https://en.wikipedia.org/wiki/Sum_of_two_squares_theorem`,
		`"Legendre's three-square theorem"`,								#Ref 3
			`https://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem`,
		`"Euler's four-square identity"`,									#Ref 4
			`https://en.wikipedia.org/wiki/Euler%27s_four-square_identity`,	
		`"Jacobi's four-square theorem"`,									#Ref 5
			`https://en.wikipedia.org/wiki/Jacobi%27s_four-square_theorem`
;
(* All numeric variables represent nonnegative integers. This is to be taken as given when reading 
 * these comments.																	*)
local
	#Efficient and cached computation of floor(sqrt(n)):
	Isqrt:= proc(n::nonnegint) option cache; local s:= isqrt(n); `if`(s^2 > n, s-1, s) end proc,

	#Cache the prime factorization of n:
	Ifactors:= proc(n::And(posint, Not(1))) option cache; ifactors(n)[2] end proc,

	#Cache representations as a sum of 2 squares:
	Sum2Sqs_cache:= proc(n::And(nonnegint, satisfies(IsSum2Sqs)))
	option cache;
		NumberTheory:-SumOfSquares(n)
	end proc,  

	#Cache the random generators:
	Rand:= proc(N::range(nonnegint)) option cache; rand(N) end proc,
	
	#Cache combinatorial counts of multinomials: 
	Multinomial:= proc() option remember; combinat:-multinomial(4, args) end proc
;
export
	#Test whether n is a sum of 2 squares:
	#n is sum of 2 squares iff every prime of the form 4*k+3 in its
	#factorization occurs to an even exponent (Ref. 2).
	IsSum2Sqs:= proc(n::nonnegint)
	option cache; 
		n > 1 implies andmap(p-> irem(p[1],4)=3 implies p[2]::even, Ifactors(n))
	end proc,

	#Express n as a sum of 2 squares, by random choice:
    Sum2Sqs:= proc(n::And(nonnegint, satisfies(IsSum2Sqs)))
    local S:= Sum2Sqs_cache(n);
        S[Rand(1..nops(S))()][] 
    end proc,
    
    #Test whether n is a sum of 3 squares: 
    #n is a sum of 3 squares iff it is NOT of the form 4^a*(8*b+7) (Ref. 3). 
    IsSum3Sqs:= proc(n::nonnegint)
    option cache;
    local r, q:= n;
        if n=0 then return true fi;
        do until irem((r:= q), 4, 'q') <> 0;
        evalb(irem(r,8) <> 7)
    end proc,

    #Express n as a sum of 3 squares, by random choice:
    Sum3Sqs:= proc(n::And(nonnegint, satisfies(IsSum3Sqs)))
    local R:= Rand(0..Isqrt(iquo(n,3))), a, b, trials;
          for trials do until IsSum2Sqs((b:= n - (a:= R())^2));
        if trials > 1 then userinfo(1, ':-Sum3Sqs', trials, "trials used for", n) fi;
          a, Sum2Sqs(b)
    end proc,

    #Euler's 4-square identity (Ref. 4):
    #This holds for any commutative ring, R. 
     #It computes C such that add(C^~2) = add(A^~2) * add(B^~2),
     #where A, B, and C are all 4-element lists.                *)
    Euler4Sq:= (
        A::[algebraic, algebraic, algebraic, algebraic],
        B::[algebraic, algebraic, algebraic, algebraic]
    )-> [
        A[1]*B[1] - A[2]*B[2] - A[3]*B[3] - A[4]*B[4],
        A[1]*B[2] + A[2]*B[1] + A[3]*B[4] - A[4]*B[3],
        A[1]*B[3] - A[2]*B[4] + A[3]*B[1] + A[4]*B[2],
        A[1]*B[4] + A[2]*B[3] - A[3]*B[2] + A[4]*B[1]
    ],

    #Express n as a sum of 4 squares, in one way (Ref. 1): 
    #Use either a straightforward algorithm or
    #distribute the work over the prime factorization and combine the results with Euler4Sq. 
    Sum4Sqs:= proc(
        n::nonnegint,
        #keyword parameter to force use of straightforward algorithm: 
        {nofactor::truefalse:= false}
    )
    local R, a, b, S, trials;
        if nofactor or isprime(n) or n < 2 then
            R:= Rand(0..Isqrt(iquo(n,4)));    
            for trials do until IsSum3Sqs((b:= n - (a:= R())^2));
            if trials > 1 then userinfo(1, ':-Sum4Sqs', trials, "trials used for", n) fi;
            S:= [a, Sum3Sqs(b)]
        else
            S:= abs~(foldr(Euler4Sq, [1,0,0,0], thisproc~((`$`@op)~(Ifactors(n)))[]))
        fi;
        sort(S)
    end proc,

    #Number of solutions, counting negatives and permutations (Ref. 5):
    #This is Jacobi's four-square theorem: The number of solutions is 8*s(n) - 32*s(n/4), where
    #        { -1/24,                    if N = 0;
    # s(N) = { sum of the divisors of N, if N::posint;
    #        { 0,                        otherwise. 
    NumSols:= proc(n::nonnegint)
    local r, q:= iquo(n,4,r);
        `if`(n=0, 1, 8*NumberTheory:-sigma(n) - `if`(r=0, 32*NumberTheory:-sigma(q), 0))
    end proc,

    #Account for multiplicity--negatives and permutations:
    #The multiplicity of a solution is 4!/mul(r!, r= sol) * 2^4/2^z, where sol counts the times each
    #number appears in a solution and z is the times 0 appears in a solution. 
    Multiplicity:= (sol::seq(nonnegint))-> 
        (T-> Multinomial(entries(T, 'nolist')) * 16/2^`if`(assigned(T[0]), T[0], 0))
            (Statistics:-Tally([sol], 'output'= table)),

    #All 4-squares solutions: 
    #Generate solutions at random until all are accounted for or a desired number are obtained.
    #If random selection produces too many duplicates, this procedure signals to Sum4Sqs to switch 
    #to the no-factor algorithm. The reason that this helps is that not every representation can be
    #found by solving for the factors and applying Euler4Sq.
    ModuleApply:= proc(
        n::nonnegint, 
        max_sols::posint:= infinity #optional argument: max # of desired solutions
    )
    local Nsols:= NumSols(n), nsols:= 0, Sols, sol, sols:= 0, duplicates:= 0, nofactor:= false;
        while nsols < Nsols and sols < max_sols do
            sol:= Sum4Sqs(n, ':-nofactor'= nofactor)[];
            if assigned(Sols[sol]) then
                ++duplicates;
                if not nofactor and (nofactor:= evalb(duplicates > sols)) then
                    userinfo(1, ':-Sums4Sqs', "Switching to no-factor algorithm at", sols, "solutions.")
                fi; 
                next 
            fi;
            Sols[sol]:= ();
            nsols+= Multiplicity(sol);
            ++sols
        od;
        userinfo(
            1, ':-Sums4Sqs', 
            sols, "solutions", `if`(nsols < Nsols, "found;", "total;"),
            duplicates, "duplicates discarded."
        );
        sort([indices(Sols)])
    end proc
;
end module:

Download Sum4Squares.mw

What you describe seems to me to be standard file buffering, not something specific to Maple. If you require access to the file from an external application (such as Windows File Explorer to see the file's size), then you can close the file after every write operation, and re-open it (in append mode) before the next write. 

Here's a much more efficient algorithm:

restart
:
#Efficient computation of floor(sqrt(n)):
Isqrt:= proc(n::nonnegint)
local s;
   for s from isqrt(n) by -1 while s^2 > n do od;
   s
end proc
:
#Test whether n is a sum of 2 squares:
IsSum2Sqs:= (n::nonnegint)-> 
   issqr(n) or andmap(p-> irem(p[1],4)=3 implies p[2]::even, ifactors(n)[2])
: 
#Express n as a sum of 2 squares:
Sum2Sqs:= proc(n::And(nonnegint, satisfies(IsSum2Sqs)))
local a:= Isqrt(n), b;
   for a from Isqrt(n) by -1 do b:= n-a^2 until issqr(b);
   a, Isqrt(b)
end proc
:
#Test whether n is a sum of three squares: 
IsSum3Sqs:= proc(n::nonnegint)
local r, q:= n;
   if IsSum2Sqs(n) then return true fi;
   do r:= q until irem(r, 4, 'q') <> 0;
   evalb(irem(r,8) <> 7)
end proc
: 
#Express n as a sum of 3 squares:
Sum3Sqs:= proc(n::And(nonnegint, satisfies(IsSum3Sqs)))
local a, b;
   for a from Isqrt(n) by -1 do b:= n-a^2 until IsSum2Sqs(b);
   a, Sum2Sqs(b)
end proc
:
#Euler's 4-square identity (true for any commutative ring): 
#Returns C such that add(C^~2) = add(A^~2) * add(B^~2) 
#where A, B, C are all 4-element lists.
Euler4Sq:= (
   A::[algebraic, algebraic, algebraic, algebraic],
   B::[algebraic, algebraic, algebraic, algebraic]
)-> [
   A[1]*B[1] - A[2]*B[2] - A[3]*B[3] - A[4]*B[4],
   A[1]*B[2] + A[2]*B[1] + A[3]*B[4] - A[4]*B[3],
   A[1]*B[3] - A[2]*B[4] + A[3]*B[1] + A[4]*B[2],
   A[1]*B[4] + A[2]*B[3] - A[3]*B[2] + A[4]*B[1]
]:
#Express n a sum of 4 squares
Sum4Sqs:= proc(n::nonnegint)
local a, b;
   if isprime(n) or n < 2 then
      for a from Isqrt(n) by -1 do b:= n-a^2 until IsSum3Sqs(b);
      return [a, Sum3Sqs(b)]
   fi;
   abs~(foldr(Euler4Sq, [1,0,0,0], thisproc~((`$`@op)~(ifactors(n)[2]))[]))
end proc
:

Example usage:

n:= rand(); Sum4Sq(n);
                       n := 927663025828
                [486168, 664104, 121392, 485318]

There is an even-more-efficient algorithm known, but the paper was behind a pay wall. (Details later).

Create your plot in Maple, for example, via

P:= plot(f, 1..3);

Then do

ExportMatrix("somefilename.csv", op([1,1], P), target= csv);

There is no (formal, legitimate) way to decrypt such files, nor should there be. One of these situations necessarily applies:

1. You want your colleagues to be able to read your source code. Then send them the source code.
2. You want your colleagues to be able to read your source code but you don't want any interlopers to be able to. Then send the source code encrypted by an external application.
3. You want your colleagues to be able to run BUT NOT read your code. Then follow Tom's Answer (send the libraries with Maple's internal encryption).

None of those scenarios require you or your colleagues to decrypt Maple's internal encryption.

You just need to add an op to your dsolve command:

dsolve({op(%), op(Initial)}, %%);

You remembered it for Initial, but not for the ODEs. The above will give you a massively complex symbolic solution, which is probably of no use to you. So, add keyword numeric and save the solution to a variable for later use:

Sol := dsolve({op(%), op(Initial)}, %%, numeric);

To address avoidance of such mistakes: I find %, %%, etc., useful for off-the-cuff computations, but I also find that they make making mistakes more likely. So, I remove them from my finished work. If the worksheet is worth saving, then it's worth deconstructing the % references (from my point of view).

Here are details regarding three points of my random prime generators---elaborate enough that I thought a new Answer is justified. The points:

1. I said before that it would be easy to modify the generators to efficiently draw multiple samples from the same interval. Actually, it's so easy that there's no point to not doing it, so both generators below operate this way.
2. I thought that it was obvious that my UniformRandPrime was uniform. Well, maybe it wasn't so obvious. Below is some emipirical evidence of uniformity.
3. It probably wasn't obvious that my procedures work quickly even when the set of primes being selected from is much too large to create. Below is evidence of that efficiency.

(Worksheet is in Reply below due to some bug in MaplePrimes.)

Your equation has two dependent (or functional) variables: x(t) and y(t). That means that there'll need to be a system of two differential equations in order to get a solution.