Beautifying is not the strength of a CAS ...

Here an old command will work:

phi[n, j] = a[j, D]-a-sqrt((-a[j, D]+a+tau[n, j])^2+4*gamma*f[c, n, j]/L[j]):
student[completesquare](%, a[j, D]);

  phi[n, j] = a[j, D] - a

           /                         2   4 gamma f[c, n, j]\1/2
         - |(a[j, D] - a - tau[n, j])  + ------------------|
           \                                    L[j]       /

PS: be aware that your symbols may have a special meaning in Maple
1) D is the differential operator, 2) gamma is Euler's constant ~ 0.57 ...

Adding (while it is not clear, over which field the task should be viewed at): in general one can not solve by radicals for degree beyond 4. There are ways using special functions. But practical it quickly becomes useless to ask for symbolic solutions. And note that even for degree=2 it is not completely trivial to get a reasonable numerical result from a symbolic solution (you will find that in "Numerical Recipies", pls search for it). Do not expect that generations chew on that and having Maple you get it solved.

hm ... I would plot. And expect a student to 'know' that cos(t) is positive for t=theta/2 between 0 and Pi/2. What a mess to use a CAS.

You may wish to care for upper and lower case in tags: matrix vs Matrix. And allow for expressions, not only single words, like cumulative distribution, numerical integration and so on

To tell the sad truth: for testing one gets 20 - 30 % less at the SW development market, documentation is even worse. And if one has it in his vita as decent task then you can imagine the rest.

Please see the help how to index a Matrix. You may read some basics (though there are some other ways to code that).

And what is 'k' in your loop?

My personal explanation is more simple: testing and documentation is boring, not repected, low paid. But not trivial at all.

http://geekandpoke.typepad.com/.a/6a00d8341d3df553ef014e8b8f1c5b970d-800wi

Perhaps it is meant to write or display erf as cumulative normal. Let me denote it by cdfN.

eval(y, erf = 't -> 2*cdfN(t*sqrt(2))-1'); 

                                                          2
             1.41421356237309 (0.0414841787338403 + 0.25 x )
        cdfN(-----------------------------------------------)
                                    x

In Maple 18 it works. You may wish to check the syntax for that command through the help (the message seems to say that one does not have the option you want to use in Maple 17)

and it shows me my answer in FF, but not in IE, still not showing the screen shot I made (visible in IE at psoting time but now it is gone). Only after logout in IE it show my sreen shot as missing.

I quit testing.

what I see is  in FF (with security settings) in IE it is ok. Have not investigated further, my knowledge on js is rather limited. But I do not allow asm javascript. You use http://www.mapleprimes.com//services/VoteService.asmx/js, no?

I will sketch it for N=2 using Digits:=15 and t for the variable.

Using series(s, t=0, 3) gives 

  -5.29235542443546*t^2+23.1633800067527*t-2.*ln(t)-102.121304227141

Since t is very small the terms at t, t^2 are certainly smaller then
DBL_EPSILON (-> google that) and adding them to the constant does not
change it numerically (it is larger than 1, DBL_EPSILON as buzz word).

Thus essentially one has to solve -N*ln(t) = eval(s, t=0), which is
given by exp. But eval(s, t=0) will throw an error (because of log),
and a correct limit would be expensive.

Now use a small number instead, t_guess = exp( eval(s, t=1e-100)/N) )
will do for a good solution. You may wish to compare it to the result
given by fsolve(s, t=0 .. infinity). Feed both into your task to check
whether results actually numerically solve the task. Allways: check it.

You should be able to refine that for your needs. That way it does not
need days, but milli seconds.

You still can proceed the save way suggested by Markiyan Hirnyk. 

NB: Neither the above nor a plain use of standard commands of a CAS
may be enough for a thesis without justifications.

Well, even as teaching assistant it would not be my job to provide that.

If someone is writing a Diploma thesis covering numerical tasks I expect that s/he either knows such things from lectures or will be able to learn about it by either asking the teaching assistant or a search engine (already Google brings up some thousand hits after entering the search term).

I will answer later.

I edited (t=l), substracting singularities is a common trick. Else proceed a Markiyan Hirnyk suggested (may be increase Digits to 15 or more) and care for difficulties beyond N = 8, may be depending on you other inputs.

Roughly: your s has a singularity -N*ln(t) and substracting it that will
be regular in t=0. As we 'know' (from the discussion) that the solution
is very close to zero and looking at the series in t=0 one can use the
value of correction:= s - singularity in t=0 to solve the equation
N*ln(t) = eval(correction, t=0), which already is good enough as solution.

Just in case the suggested way is not fast enough ... Also note that the
problem can be solved up to N=8. Beyond that I have doubts that a real
solution exists.

Edited & PS: 
I used t instead of l (the letter, not the number 1 or sqrt(-1) = I) for readability