To get both symbols and connecting lines, I think you need to use separate pointplot commands (or a plot and a pointplot), and then you can combine them using display.  Thus:
 

> display([pointplot(stressStrainData, connect=true, colour=red), 
           pointplot(stressStrainData,  symbol=circle, colour=blue)], labels=[engrStrain, engrStress]);

It is true that assigning to a variable is bad if you're using a subscripted version of that variable.  But that's not what's happening in your original example, or my simplified version

> simplify(sqrt(omega^2 + omega[0]));

Error, (in content/content) invalid arguments

This is a bug.  It shouldn't happen.  If Maple allows you to create the expression sqrt(omega^2 + omega[0]), then simplify should be able to operate on it without producing a cryptic error message.  The only situations where simplify should result in an error is when there is something wrong with the expression (e.g. a fraction where the denominator simplifies to 0).

It is true that assigning to a variable is bad if you're using a subscripted version of that variable.  But that's not what's happening in your original example, or my simplified version

> simplify(sqrt(omega^2 + omega[0]));

Error, (in content/content) invalid arguments

This is a bug.  It shouldn't happen.  If Maple allows you to create the expression sqrt(omega^2 + omega[0]), then simplify should be able to operate on it without producing a cryptic error message.  The only situations where simplify should result in an error is when there is something wrong with the expression (e.g. a fraction where the denominator simplifies to 0).

Something like this, perhaps?

with(LinearAlgebra): 
Totals:= [0,0,0,0]:
to 100 do
   A := RandomMatrix( 2, 2, generator=-100..100 );
   E := Eigenvalues( A );
   if   has(E,I)              then Totals[1]:= Totals[1]+1
   elif has(E,0)              then Totals[2]:= Totals[2]+1
   elif simplify(E[1]*E[2])<0 then Totals[3]:= Totals[3]+1
   elif simplify(E[1]*E[2])>0 then Totals[4]:= Totals[4]+1
   else print("uh-oh, something's wrong")
   end if;
 end do:
 Totals/100.0;

Something like this, perhaps?

with(LinearAlgebra): 
Totals:= [0,0,0,0]:
to 100 do
   A := RandomMatrix( 2, 2, generator=-100..100 );
   E := Eigenvalues( A );
   if   has(E,I)              then Totals[1]:= Totals[1]+1
   elif has(E,0)              then Totals[2]:= Totals[2]+1
   elif simplify(E[1]*E[2])<0 then Totals[3]:= Totals[3]+1
   elif simplify(E[1]*E[2])>0 then Totals[4]:= Totals[4]+1
   else print("uh-oh, something's wrong")
   end if;
 end do:
 Totals/100.0;

Assuming you're talking about the differential equation system x' = A x, I think you want to change that "isolated" to "non-isolated".
 

Assuming you're talking about the differential equation system x' = A x, I think you want to change that "isolated" to "non-isolated".
 

Or you might try

> ArrayTools[Reshape](Statistics[Sample](Normal(0,1),4),2,2);

Or you might try

> ArrayTools[Reshape](Statistics[Sample](Normal(0,1),4),2,2);

The normal way to define a function is using  ->.  With N40(xi):= ...  you are literally defining the value of the procedure N40 only on the literal input xi (technically, you are assigning an entry in the remember table of N40), and this will not affect N40(anything else).
 

If you are using 2D input, N40(xi):= ... will cause a popup window asking you to clarify what you mean, a function definition or a remember table assignment.  If you choose function definition, the result will be as if you used ->.  But you might as well use -> in the first place.

The normal way to define a function is using  ->.  With N40(xi):= ...  you are literally defining the value of the procedure N40 only on the literal input xi (technically, you are assigning an entry in the remember table of N40), and this will not affect N40(anything else).
 

If you are using 2D input, N40(xi):= ... will cause a popup window asking you to clarify what you mean, a function definition or a remember table assignment.  If you choose function definition, the result will be as if you used ->.  But you might as well use -> in the first place.

This is a different problem.  You made c a local variable of NCcoef.  The c[0], ..., c[3] returned by your first NCcoef are referring to this local variable.  In your first A, on the other hand, you are using global c[0],..., c[3].  Although the names look the same to you, Maple considers them to be completely different, so eval fails to do any substitution.  It would work if you removed c from the local variables declaration.

Your second version succeeds because names constructed using || are always global, not local. 

In Maple 13 I get:

> dsolve({ics, ode});

theta(t) = -2*arctan(12*(exp(7/25*t+2*I*Pi*_Z3-4*I*Pi*_Z4)-1)/(-16+9*exp(7/25*t+2*I*Pi*_Z3-4*I*Pi*_Z4)))

> simplify(%);

theta(t) = -2*arctan(12*(exp(7/25*t)-1)/(-16+9*exp(7/25*t)))

Note however that this solution is only valid until t = 25/7*ln(16/9) = 2.054871946.  A better solution would use the two-argument form of arctan:

theta(t) = -2*arctan(-12*exp(7/25*t)+12,16-9*exp(7/25*t)) 

In Maple 13 I get:

> dsolve({ics, ode});

theta(t) = -2*arctan(12*(exp(7/25*t+2*I*Pi*_Z3-4*I*Pi*_Z4)-1)/(-16+9*exp(7/25*t+2*I*Pi*_Z3-4*I*Pi*_Z4)))

> simplify(%);

theta(t) = -2*arctan(12*(exp(7/25*t)-1)/(-16+9*exp(7/25*t)))

Note however that this solution is only valid until t = 25/7*ln(16/9) = 2.054871946.  A better solution would use the two-argument form of arctan:

theta(t) = -2*arctan(-12*exp(7/25*t)+12,16-9*exp(7/25*t)) 

> is(0< mu<1);

(entered in 2D Math input) produces the result false  (i.e. Maple is just telling you that this does not follow from the assumptions you've made).  I think you want

> assume(0 < mu < 1);

But as for your error: it looks like you've assigned a value to theta(t)  in the line

theta(t):= ...

Once you've done that, you don't have theta(t) as an unknown function for your ODE any more.