The fisher information matrix is:
_____________________________
I(u,s) =  1/n(3s^2)  0 

 0  (pi^2 + 3) / n(9s^2) 

Where n is the number of sample points.
Variance covariance matrix is the inverse of the Fisher Information Matrix... I^1
You can then use that to calculate the confidence intervals on u and s using the positions on the I^1 matrix
Var(u) = I^1(1,1) %Row 1 Col 1
Var(s) = I^1(2,2) %Row 2 Col 2
Where s is the point estimate parameter from MLE.
Where p is your confidence number (eg. 0.9)
Then you can assume the confidence intervals follow a normal distribution (not always true unless you do inf samples but a quick and dirty way to do it):
u_low = u  InvNorm[ (1+p)/2] sqrt[ Var(u)]
u_upper = u + InvNorm[ (1+p)/2] sqrt[ Var(u)]
For s (noting that it must be greater then 0) it can be assumed that this is a truncated normal distribution:
s_lower = s. exp{ InvNorm[ (1+p)/2] sqrt[ Var(u)] / s }
s_upper = s. exp{ InvNorm[ (1+p)/2] sqrt[ Var(u)] / s } note not negative on denominator
From there you have the confidence intervals of your parameter estimates.
Now you can do the pdf, cdf etc confidence intervals using a very quick monte carlo at each point.
I'm sure there is a quicker way though...
Good luck..
Occa
