You have to extract the z-values of the plot, and calculate the corresponding x- and y-values.
I give an example

restart;
p := plot3d(x^2-y^2, x=-1..1,y=-1..2):
pdata := plottools:-getdata(p); zvals := pdata[3]:
# the x- and y-range are divided in 48 equal intervals.
# we compute the x- and y-values in this grid:
h := 2/48: xvals := seq(evalf(h*i),i=-24..24);
h := 3/48: yvals := seq(evalf(h*i), i=-16..32);
# an arbitrary point on the graph:
[xvals[9],yvals[23],zvals[9,23]];
# check:
xvals[9]^2-yvals[23]^2;

I have put your data in a textfile: HassanData.txt
Now you can read this in a n×4-Matrix, and construct two new n×2-Matrices consisting of the x-values as the first column, and  as the second column the difference of the exact value and respectively the CN and the CS-value:

A := Matrix(readdata("HassanData.txt",4));
CN := <A[..,1]|A[..,4]-A[..,2]>;
CS := <A[..,1]|A[..,4]-A[..,3]>;
plot( [CN,CS] );

Maple doesn't recognize sqrt(2) as a positive number.
The best you can do is to copy the actual arguments to local variables and test if these are all positive reals.
Due to rounding errors the case (a+b)=c must be replaced by an inequality condition.
So you get

isTriangle:=proc(a,b,c)
  local A,B,C;
  A,B,C := seq( evalf(args[i]), i=1..3 );
  if not(type(A,positive) and type(B,positive) and type(C,positive)) then 
    error "not all arguments are positive constants"
  end if;
  if (A+B) > C then print (a,b,c,`are lengths of a triangle`)
#this is according to the triangle inequality theorem
  elif (A+B) < C then print (a,b,c,`are not lengths of a triangle`);
  elif abs(A+B-C) < 1e-10 then print (a,b,c, `does not adhere to the triangle inequality theorem`)
  end if
end proc;

I take the constant 1 (instead of 10^10). Then I get:

DE := diff(z(t),t,t) = -((sqrt(diff(z(t),t)^2+z(t)^2)*diff(z(t),t))+z(t));
DE2 := diff(z(t),t)=y(t),  subs( {diff(z(t),t)=y(t)}, DE );
with(DEtools):
DEplot( {DE2}, [z,y], t=0..10, 
   {[z(0)=-2,y(0)=2],[z(0)=-1,y(0)=-2],[z(0)=2,y(0)=-2]}, 
   z=-2..2, y=-2..2, arrows=smalltwo, linecolor=blue );

which looks like figure 3.5

Your formula is still not unambiguous (unbalanced parentheses). I suppose that you mean:

DE := diff(z(t),t,t) = -((sqrt(diff(z(t)^2,t))+10^10*diff(z(t),t))+10^10*z(t))/diff(z(t),t);

Now rewrite to a system of two first order equations:

 DE2 := diff(z(t),t)=y(t),  subs( {diff(z(t),t)=y(t)}, DE );

A phase portrait with two trajectories:

with(DEtools):
DEplot( {DE2}, [z,y], t=0..1, 
   {[z(0)=5e9,y(0)=2e10],[z(0)=5e9,y(0)=1.5e10]}, 
   z=5e9..2e10, y=5e9..2e10, arrows=smalltwo, linecolor=blue );

By the way, because of the huge constants in the equation, you can omit the sqrt term and simplify to

DE3 := diff(z(t), t) = y(t), diff(y(t), t) = -(10^10*y(t)+10^10*z(t))/y(t);

which gives you the same phase portrait.

If you can live with 1.10^-8:

textplot( [1,1,typeset(beta=1e-8)] );

Or, without the 1. before the 10^-8

textplot( [1,1,typeset(beta=` 10`^` -8`)] );

(mind the spaces!)

I think that my first answer wasn't clear enough.

So this procedure prints two plots (a print results in no "output"!) and outputs only the number 56.
You have to change something in the procedure: assign the plot to an output variable -- see my first answer!
See also: ?paramater modifiers

Download rk2b_sys_Task34.mw

You have to make the plot the (or an) output of the procedure: so the plot command should be the last statement before end proc.
But perhaps the procedure has other relevant output. Then you can assign the plot to an output variable:

f := proc( input parameters, P::evaln )
  ...
  P := plot(...)
  ...
  ...
end proc;


f(inputs, FirstPlot );

plots:-display( {FirstPlot, other plots} );

Your g[7] is not calculated.
The innerproduct is an easy integral, so you don't need evalf(Int(...
You can do evalf when you calculate the g[k]'s.

This is a revised worksheet:

gramschmidt2_1A.mw

This equation can be solved exactly:

restart;
cons := {c1=(-21+.7*I)/2.3, c2=(2.3*6.494)*10^13, c3=(-21+.7*I)*6.494*10^13 }:
m := 1 + c1*sqrt(c2-x^2)/sqrt(c3-x^2):
s := solve(m,{x}):
X := eval( s[1], cons );
           
            { x = 1.295020172*10^7  + 26513.56567*I }
           
This is indeed inside your rectangle.
However:

radsimp(eval(m,s[1]));

                               2

Wrong answer!
But when we rewrite the equation for m, we get

 m := 1 + c1*sqrt((c2-x^2)/(c3-x^2)):
radsimp(eval(m,s[1]));
                               2
# also false.
  ...and pulling a minus sign outside the sqrt:
m := 1 + c1*I*sqrt((-c2+x^2)/(c3-x^2)):
radsimp(eval(m,s[1]));
                               0

This strange behaviour has to do with the branch that is chosen bij the sqrt function: sqrt(a/b) is not always equal to sqrt(a)/sqrt(b):

 a,b := exp(2/3*Pi*I),exp(-2/3*Pi*I):
evalc(sqrt(a/b)), evalc(sqrt(a)/sqrt(b));
                1   1    (1/2)    1   1    (1/2)
                - - - I 3     , - - + - I 3     
                2   2             2   2         

M := n -> Matrix( [seq( [f(i),g(i),h(i)], i=1..n )] ):
M(4);

These are a system of nonautonomous ODE's. I try to find a numeric solution of the IVP:

t := diff(X(x), x) = -(1-6*R(x)^(1/2))^(1/2)*x*X(x)/(X(x)*x+R(x))^(1/2), 
     diff(R(x), x) = (1-6*R(x)^(1/2))^(1/2)*x^2*X(x)/(X(x)*x+R(x))^(1/2):
init := R(0) = 0, X(0) = 100:
sol := dsolve( {t,init}, {R(x),X(x)}, type=numeric, range=0..2 ):
Warning, cannot evaluate the solution past the initial point, problem may be complex, 
initially singular or improperly set up

Indeed, "initially singular": both denominators are zero at x=0.
So I take other IV's:

init := R(1) = 0, X(1) = 100:
sol := dsolve( {t,init}, {R(x),X(x)}, type=numeric, range=1..1.05 ):
plots:-odeplot(sol, [x,R(x)], refine=50);

plots:-odeplot(sol, [x,X(x)], refine=50);

(Edited 1/29/2013): You want the units of the x axis to be x/M:

plot( int(t^2*exp(Pi*t/10^12*M), t=0..x ), x = 0 .. 1000, 
  tickmarks=[[seq(100*i=.1*i,i=1..10)],default], 
  labels=[typeset(` x/M`*` 10`^` -5`),""] );

Try to get your tables as a list of lists: [[t1,x1,z11],[x2,t1,z21],...]
Use ?surfdata.

simplify( f, {A^2+1=M^2} ) assuming M>0;
                              (1/2)     
                          / 2\         2
                       -2 \M /      + M 

Strange, that another simplification is needed for sqrt(M^2):

simplify( simplify( f, {A^2+1=M^2} ) ) assuming M>0;
                                   2
                           -2 M + M