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.

The line with the add commands should be (at least in the current version of Maple):

L:= (1/3)*h*(y(x[1])+4*add(y(x[i]), i = 2 .. N, 2)+2*add(y(x[i]), i = 3 .. N-1, 2)+y(x[N+1]));

The result, after evalf, should be 4354.742601 (as you got). If you get something else (not close) or an error, use this:

L:= (1/3)*h*(y(x[1])+4*add(y(x[2*i]), i= 1 .. N/2)+2*add(y(x[2*i+1]), i= 1 .. N/2-1)+y(x[N+1]));

I'm sure that this latter syntax will work in Maple 16. (Current Maple allows a third argument for add--the "skip"--just as you used h as a third argument for seq.)

You may compare your results with:

Student:-Calculus1:-ApproximateInt(y(X), X= a..b, method= simpson, partition= N);
evalf(%);

Also note that your function y(X) is equivalent to a hyperbolic cosine cosh(X/C).

Unit should be capitalized, not unit. But you probably don't need to set anything to access all 32 G from Maple.

Write each method as a procedure (proc(...)). Let's say that these procedures are named Jacobi, GaussSeidel, and SOR. For each, do 

CodeTools:-Usage(Jacobi(...)),

etc.

This returns the result of the inner procedure. The times and memory are printed as side effects.

If you're timing things that take "minutes", this should be sufficiently accurate. Some special care is needed to time very small things (say, < 0.2 seconds).

Algorithm:

1. Select points Xj at random in C. Discard any at distance <= R1 from X1. Continue until you have m-1 points.
2. Compute min pairwise distance among X2, ..., Xm. Call this eps1.
3. Compute min({|X1-Xj| - R1 : 2 <= j <= m}). Call this eps2.
4. Let all of R2, ..., Rm = min(eps1, eps2) / 3.

For each instance of beta(10) you need to put a space between the beta and the left parenthesis. 

The appropriate command is add, not sum. It seems that everyone wants to use sum, which is primarily for infinite, symbolic, and/or indefinite summation. If you have specific integer values of the index, use add. Yes, you could use quote marks, as suggested by NM, but the better solution is to use add.

Also, you should use eval instead of subs. It seems that everyone wants to use subs, a command which has no knowledge of mathematics.

add(eval(diff(x^2, [x$k]), x= 0), k= 1..2);